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Abstract 


This  report  is  the  second  in  a  series  of  reports  discussing  mitigation  of  blasts  using  water  mists. 
The  previous  report  described  numerical  simulations  of  a  TNT  blast.  The  present  report  has  combined 
that  numerical  procedure  with  a  two-continuum  model  for  dispersed-phase  calculations.  The  specific 
approach  used  for  the  dispersed-phase  calculations  is  the  sectional  approach,  and  is  described  in  this 
report  along  with  a  series  of  simulations  conducted  for  validation,  and  also  for  the  purpose  of  better 
understanding  interactions  between  shock  waves  and  particle-seeded  gases. 

A  new  derivation  of  the  sectional  approach  is  presented  in  this  report  to  provide  details  and  to  outline 
the  limitations  of  this  approach.  This  approach  is  then  incorporated  into  a  flux-corrected-transport  high 
speed  flow  code,  and  is  used  to  examine  mitigation  characteristics  of  glass  particles  and  water  droplets  on 
impinging  shock  waves.  First  a  series  of  validation  cases  are  presented  colnparing  simulations  against 
available  experiments.  Results  show  that  the  present  solution  technique  accurately  predicts  the  shock 
wave  behavior  in  particle-seeded  gases.  Further  simulations  were  conducted  to  investigate  the  effect  of 
particle  size  and  vaporization  for  both  non-decaying  and  decaying  shocks.  For  steady  shock  waves,  it 
was  determined  that  for  finite-sized  particles,  the  shock  front  pressure  and  maximum  overpressure  are 
not  equal,  with  the  shock  front  pressure  typically  being  smaller  than  the  unmitigated  shock  pressure,  and 
the  maximum  overpressure  being  greater  than  the  unmitigated  shock  pressure  for  non-decaying  shocks, 
dependent  primarily  on  the  mass-loading  of  the  particles.  For  decaying  shock  waves,  there  is  an  optimum 
particle  size  which  gives  the  most  effective  mitigation  for  the  peak  overpressure,  but  is  specific  to  the 
shock  tube  geometry  considered.  Vaporization  effects  tended  to  be  minimal  for  the  cases  presented  in 
this  report,  because  the  shock  waves  were  fairly  weak  resulting  in  only  moderate  temperatures  for  the 
shocked  driven  gases.  Blasts  tend  to  be  similar  to  the  decaying  shock  wave  solutions,  except  that  they 
are  considerably  more  powerful,  and  should  have  high  enough  temperatures  such  that  vaporization  will 
play  a  more  significant  role  in  mitigation. 
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BLAST  MITIGATION  BY  WATER  MIST 
2)  SHOCK  WAVE  MITIGATION  USING  GLASS  PARTICLES 
AND  WATER  DROPLETS  IN  SHOCK  TUBES 

1  Introduction 

Explosions  and  fires  have  always  represented  a  danger  to  naval  ships.  Large  amounts  of  fuel  and  ord¬ 
nance  on-board  these  ships  provide  plentiful  fuel  for  explosions  and  fires  even  in  peacetime.  Recent  events 
such  as  the  USS  Cole  have  highlighted  the  need  to  be  able  to  mitigate  the  effects  of  explosions  on  board 
ships,  whether  caused  by  internal  or  external  factors. 

Several  recent  studies  have  looked  at  water  to  mitigate  the  effects  of  explosions.  Using  water  is  attractive 
for  several  reasons.  The  large  heat  capacity  and  latent  heat  of  vaporization  allows  water  to  easily  absorb 
a  large  amount  of  energy.  In  addition,  water  systems  have  a  great  deal  of  flexibility  and  low  cost,  are 
environmentally  safe  and  clean  to  use,  and  can  serve  a  dual  role  of  fire  suppression  in  peacetime  situations 
and  explosion  mitigation  during  war  and  terrorist  situations.  A  recent  summary  of  efforts  to  use  water  and 
water  mist  systems  to  mitigate  the  effects  of  blasts  and  explosions  is  given  in  [1]. 

Of  particular  interest  for  the  current  project  is  mitigation  of  the  initial  and  quasi-steady  overpressures  felt 
within  an  enclosure  that  is  subjected  to  an  explosive  blast.  Previous  studies  [2]  demonstrated  that  the  blast 
characteristics  within  simple  enclosures  for  2.12  kg  of  TNT  could  be  effectively  simulated,  reproducing 
overpressures  and  quasi-static  pressures  that  were  similar  to  experimental  results  [3].  The  next  step  is 
to  develop  a  spray  model  and  incorporate  it  into  the  blast  simulation,  to  determine  the  amount  of  actual 
mitigation  effected  by  the  water  spray. 

There  has  been  a  substantial  amount  of  research  on  spray  models  within  the  literature.  Sirignano  [4]  pro¬ 
vides  a  review  of  currently  popular  models  and  details  of  useful  solution  procedures.  However,  the  majority 
of  these  spray  and  particle  models  have  been  developed  for  low  speed  flows  such  as  occurring  in  engines 
or  fire  situations.  There  is  unfortunately  relatively  little  information  on  the  effect  of  an  impinging  shock 
wave  on  vaporization  of  water  or  fuel  droplets,  although  some  work  has  been  done  investigating  droplet 
breakup  processes  when  exposed  to  an  impinging  shock  wave  [5,  6].  These  results  are  for  substantially 
larger  droplets  than  our  current  interest. 

The  intent  of  this  report  is  to  describe  the  spray  model  that  will  be  used  for  the  subsequent  blast  mitiga¬ 
tion  studies,  and  use  this  model  to  characterize  shock  tube  flow  where  the  driven  section  is  seeded  with  glass 
particles  and  water  droplets.  One-dimensional  shock  tube  simulations  were  chosen  because  shock  tubes  have 
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many  similarities  with  explosive  blasts,  but  are  simple  enough  that  concurrent  experimental  studies  can  also 
be  done  by  interested  researchers.  In  addition,  some  experimental  studies  with  particles  have  already  been 
done  without  vaporization  [7].  Of  particular  interest  for  our  numerical  simulations  is  the  reduction  in  shock 
wave  speed  and  the  maximum  overpressure  felt  for  different  configurations  at  different  locations  in  the  shock 
tube.  These  are  the  two  main  quantities  that  can  be  easily  compared  with  experimental  results. 

2  Transmission  of  Shock  Waves  through  Particle-Seeded  Gases 

Due  to  the  importance  of  “dusty”  gases  (gases  seeded  with  particles  that  occupy  only  a  small  fraction  of 
the  volume  but  have  a  significant  mass),  there  has  been  a  considerable  amount  of  research  on  the  transmis¬ 
sion  of  shock  waves  through  particle-seeded  flows.  For  this  report,  we  are  interested  in  how  particles  and 
droplets  affect  the  pressure  and  temperature  developed  behind  the  shock  wave  as  well  as  the  velocity  of  the 
shock  wave. 

The  first  step  is  to  understand  how  the  particles  affect  a  typical  shock  wave  as  it  penetrates  from  an 
unseeded  (pure)  gas  into  a  particle-seeded  gas.  The  shock  wave  essentially  hits  a  material  interface  as  it 
passes  into  the  particle-seeded  gas.  As  happens  at  most  material  interfaces,  depending  on  the  properties  of 
the  particles,  the  shock  wave  is  transmitted  through  the  interface  and  a  reflected  wave  propagates  back  into 
the  initial  (unseeded)  mixture  as  shown  in  Figure  1  for  infinitely  small  particles.  Depending  on  the  material 
properties  of  the  gas  mixtures  on  either  side  of  the  reflected  wave,  the  wave  can  either  be  a  compressive 
wave  {U2>  <  U2)  or  a  rarefaction  wave  (u2'  >  U2).  For  the  types  of  particles  considered  in  this  paper, 
compressive  reflected  shocks  are  generated  resulting  in  lower  gas  velocities  U2'  <  U2  and  higher  pressures 
P21  >  P2-  For  particles  of  finite  size,  the  a:  —  <  diagram  shown  in  Figure  1  requires  some  slight  modifications 
that  are  discussed  later.  We  are  interested  in  determining  the  pressme  P2'  between  the  reflected  wave  and 
the  transmitted  shock  wave. 

For  particles  that  do  not  evaporate  or  sublimate,  a  useful  limit  that  allows  analysis  of  the  situation 
shown  in  Figure  1  is  the  limit  as  the  particles  become  infinitely  small.  In  this  case,  the  particles  remain  in 
equilibrium  with  the  surrounding  gas  and  an  “equivalent”  gas  formulation  can  be  used  [8].  For  the  equivalent 
gas,  the  total  density  is  represented  as  pt  =  p(l  -t-  rj),  where  p  is  the  gas  density,  and  7/  is  the  mass-loading. 
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Figure  \\  x  —  t  diagram  of  shock  wave  from  an  unseeded  (pure)  gas  impinging  on  a  particle-seeded  gas. 
The  dashed  line  shows  the  material  interface  separating  the  unseeded  and  the  seeded  gas,  and  the  dotted  line 
is  a  continuation  of  the  incident  shock  wave  for  the  unseeded  gas. 

The  pressure  is  computed  using  the  ideal  gas  law  with  the  gas  density  and  temperature,  and  the  enthalpy  of 
the  equivalent  gas  is  given  hy  H  =  p{\  +  r}c)CpT,  where  Cp  is  the  gas  specific  heat,  and  c  is  the  ratio  of 
specific  heats  between  the  solid  particles  and  gas  (c  =  C^/Cp). 

By  using  these  definitions,  we  develop  expressions  for  the  (P,  v)  Hugoniot  and  the  (P,  u)  Hugoniot  for 
particle-seeded  flows  that  are  useful  for  understanding  the  transmission  of  shocks  through  material  inter¬ 
faces,  where  v  is  the  specific  gas  volume  {v  =  1/p).  The  Hugoniot  expressions  are  developed  assuming 
a  single  steady  shock  wave  and  considering  mass,  momentum,  and  energy  conservation  through  the  shock. 
For  ideal,  calorically  perfect  gases  (meaning  that  the  specific  heat  is  constant  and  independent  of  temperature 
and  pressure),  the  (P,  v)  Hugoniot  for  a  particle-seeded  flow  is: 

-I-  ?7c)[P2W2  “  Pl^l]  =  ^(P2  “  Pl){vi  +  V2)  (1) 

where  7  is  the  specific  heat  ratio  of  the  pure  gas.  Pi,  ui  are  the  unshocked  pressure  and  gas-phase  specific 
volume  respectively,  and  P2,  V2  are  the  shocked  pressure  and  gas-phase  specific  volumes.  Another  Hugoniot 
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expression  useful  for  shock  wave  interactions  with  material  interfaces  is  the  (P,  u)  Hugoniot,  where  u  is  the 
velocity  in  the  lab  reference  frame: 


1  +  ^,  _  12 

{A  +  l/2){P2-Pi)+APi  >ltii  '  ^  ® 

where  A  =  :^(1  +  rjc).  These  quantities  can  be  related  to  the  shock  wave  velocity  through  the  Raleigh 
line,  which  is  developed  from  a  manipulation  of  the  mass  and  momentum  conservation  equations: 


7M2 


1  \P2/Pi-1 
1+riJ  1  -  v^/vi 


(3) 


where  M  is  the  Mach  number  of  the  shock  wave  with  respect  to  the  pure  unshocked  gas.  The  Hugoniot 
expressions  along  with  the  Raleigh  expression  are  extremely  useful  when  determining  the  characteristics  of 
shocks.  Using  these  expressions,  one  can  determine  the  pressure  increase  or  decrease  resulting  from  a  shock 
wave  impinging  on  a  particle-seeded  gas,  knowing  that  the  pressure  and  velocity  in  section  2’  and  3  from 
Figure  1  are  equal. 

Results  from  this  analysis  are  shown  in  Figure  2.  This  figure  specifically  looks  at  the  effect  of  seeding 
with  glass  particles  in  air.  Notable  among  the  results  is  that  the  shock  pressure  is  actually  increased  for 
the  seeded  flow  (Fig.  2b)  while  the  velocity  of  the  shock  wave  and  the  shocked  gases  (Figs.  2a  and  2d) 
decreases.  Non-vaporizing  water  droplets  have  a  similar  impact  on  the  shock  wave  velocity  and  shock 
pressure;  however,  the  effect  of  vaporization  is  difficult  to  account  for  in  such  a  simplified  analysis. 

Particle  seeding  of  gases  to  mitigate  the  pressure  associated  with  steady  shocks  is  shown  to  be  ineffectual 
from  the  above  plots.  However,  explosions  are  typically  more  complex  than  the  simple  steady  shock.  A 
representative  geometry  that  more  closely  approximates  an  actual  blast  is  a  shock  tube  experimental  setup 
with  a  small  driver  section.  Shock  tubes  have  been  used  extensively  to  understand  different  aspects  of 
shock  waves  and  rarefaction  waves.  They  provide  an  ideal  experimental  setup  because  the  flow  is  largely 
one-dimensional  and  thus  facilitates  mathematical  modeling. 


A  shock  tube  is  made  up  of  a  cylindrical  or  rectangular  pipe  divided  into  two  compartments  separated  by 
a  diaphragm  of  thin  material  as  shown  in  Figure  3.  The  smaller  compartment  contains  gas  at  high  pressure 
and  is  called  the  driver  section.  The  larger  compartment  is  called  the  driven  section  and  contains  a  gas  at 
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(C)  (d) 

Figure  2:  Mach  number  (a),  pressure  (b),  density  (c),  and  gas  velocity  ratios  (d)  are  given  for  a  shock 
wave  impinging  on  a  seeded  gas  as  a  function  of  the  shock  wave  Mach  number  in  unseeded  gas.  Assumes 
infinitely  small  glass  particles  with  Cs  =  766  J/kg  K. 
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:  Driver  section 

Seeded  driven  section 

!Pi=4.052bar,T=li92K 

k 

,  P=1.013bar,T=:298K 

Unseeded  driven  section 
P=:1.013bar,  T=298  K 


Figure  3:  Schematic  of  experimental  setup  of  Sommerfeld  [7]. 

lower  pressures.  The  device  is  operated  by  bursting  the  diaphragm,  at  which  point  the  driver  gas  expands 
out  (much  like  a  piston),  creating  a  shock  wave  that  propagates  through  the  driven  gas  while  a  rarefaction 
wave  travels  through  the  driver  gas.  The  shock  tube  experiments  that  we  are  interested  in  further  divide  the 
driven  section  into  an  unseeded  region  and  a  seeded  region.  Of  the  shock  tube  experiments  with  particles 
and  droplets,  the  Sommerfeld  experiments  [7]  remain  the  most  useful  for  explaining  the  various  effects  of 
the  particle  on  the  shock  structure. 

Figure  4  shows  sax  -  t  diagram  of  a  shock  tube  experiment  with  a  seeded  region  in  the  driven  section, 
as  shown  in  Figure  3.  The  first  event  that  occurs  is  the  bursting  of  the  diaphragm  at  f  =  0.  This  generates 
a  shock  wave  that  propagates  through  the  driven  gases  (So),  and  a  rarefaction  wave  that  propagates  through 
the  driver  gases  (R).  The  material  interface  between  the  driver  and  driven  gases  moves  forward  at  a  slower 
speed  and  is  labeled  C.  Eventually  the  front  shock  wave  impinges  on  the  particle-seeded  gases.  Again, 
this  results  in  a  transmitted  shock  wave  (S)  and  a  reflected  wave  (Sr).  The  interface  between  seeded  and 
unseeded  driven  gas  is  labeled  P.  The  slope  of  the  material  interfaces  P  and  C  are  equal,  indicating  both 
move  at  the  same  velocity  (u2>  in  Figure  1). 

For  finite  sized  particles,  a  relaxation  period  occurs  before  the  shock  wave  reaches  an  equilibrium  struc¬ 
ture  where  it  again  propagates  through  the  medium  at  a  constant  velocity.  This  is  different  than  the  particle- 
gas  equilibrium  mentioned  previously,  because  when  the  shock  structure  reaches  equilibrium,  not  all  the 
particles  within  the  gas  are  necessarily  in  equilibrium  with  their  surrounding  environment.  It  is  an  equilib¬ 
rium  in  the  sense  that  the  shock  stmcture  does  not  change  unless  it  interacts  with  an  outside  force.  The  region 
after  the  shock  wave  impinges  on  the  seeded  flow  but  before  the  equilibrium  shock  structure  is  reached  is 
called  the  transition  region  of  the  shock.  The  effect  of  this  relaxation  period  on  the  location  of  the  seeded- 
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Figure  4:  x  —  t  diagram  of  a  shock  tube  experiment.  The  So  line  indicates  the  initial  incident  shock  wave, 
S  the  transmitted  shock,  Sr  is  the  reflected  compressive  waves,  Jf2  is  the  expansion  fan  initiated  at  the 
diaphragm,  Co  and  C  is  the  contact  discontinuity  separating  driver  and  driven  gases,  and  P  is  the  material 
interface  separating  unseeded  and  seeded  driven  gases.  The  solid  lines  represent  mitigation  with  infinitely- 
small  particles,  the  dashed  lines  represent  mitigation  with  finite-sized  particles,  and  the  dotted  lines  represent 
the  unmitigated  solution. 

unseeded  interface  (P)  and  the  front  shock  wave  (5)  are  shown  as  dashed  lines  in  Figure  4.  Note  that  after 
this  transition  period,  the  shock  wave  velocity  and  seeded-unseeded  interface  velocity  for  finite  particles 
approach  the  same  velocities  as  for  infinitely  small  particles,  however,  the  material  interface  is  displaced 
closer  to  the  initial  contact  discontinuity. 

This  transition  region  is  important  to  understand  because  it  controls  how  far  the  particles  penetrate 
into  the  unseeded-driven  gases  and  whether  they  penetrate  into  the  driver  gases.  This  can  considerably 
influence  the  mitigation  aspects  for  droplets,  because  penetration  into  the  high  temperature  driver  gases  will 
result  in  considerably  more  vaporization  and  reduce  the  force  with  which  the  driver  gases  are  pushing  the 
driven  gases.  Therefore  accurately  modeling  this  transition  region  is  essential  to  give  us  confidence  in  our 
simulations. 

In  shock  tubes  with  short  driver  sections  and  long  driven  sections,  the  rarefaction  wave  reflects  off  of 
the  end  of  the  tube,  and  then  propagates  toward  the  initial  shock.  When  it  arrives  at  the  shock,  it  relieves 
the  shock  pressure,  slowing  down  the  shock  wave  and  eventually  causing  the  shock  wave  to  dissipate.  This 
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situation  more  closely  mirrors  an  explosive  situation,  where  the  “driver”  region  (the  volume  of  the  original 
explosive)  is  very  small  and  the  rarefaction  wave  almost  immediately  begins  to  relieve  the  pressure  at  the 
shock  front.  For  that  reason,  we  are  interested  in  examining  cases  both  with  long  driver  sections  and  short 
driver  sections. 

For  the  simulations  conducted  within  this  report,  the  shock  tube  is  set  up  similar  to  the  Sommerfeld 
experiments  as  is  shown  in  Figure  3.  For  the  experiments,  the  driver  section  was  2  m  long,  the  unseeded 
driven  section  was  1.05  m  long,  and  the  seeded  driven  section  was  4.76  m  long.  The  seeded  region  is 
uniformly  seeded  with  glass  particles  at  a  specified  mass-loading.  Similar  geometries  are  also  introduced  to 
study  different  aspects  of  particle-  and  droplet-shock  wave  interactions  that  are  not  as  clear  with  the  original 
Sommerfeld  geometry. 

3  Model  Development 

The  approach  chosen  for  modeling  sprays  is  the  sectional  Eulerian  approach  introduced  by  Tambour  [9], 
and  expanded  on  by  Tambour  [10, 11]  and  other  researchers  [12, 13].  Our  research  group  has  considerable 
experience  with  using  this  approach  for  low  speed  fire  mitigation  studies  [14-17] ,  where  it  has  worked  quite 
well.  Our  purpose  for  this  paper  is  to  extend  this  approach  for  high  speed  flows  such  as  found  in  shock  tubes 
and  blast  simulations. 

Lagrangian  tracking  of  particles  and  droplets  has  been  widely  used  to  simulate  these  types  of  flows 
because  of  its  flexibility  and  accuracy.  However,  to  do  Lagrangian  tracking  for  individual  particles  is  very 
computationally  very  prohibitive  for  large  numbers  of  droplets.  The  Eulerian  sectional  approach  provides  an 
alternate  method  for  doing  these  computations  that  is  well  suited  to  these  problems.  There  are  two  general 
assumptions  that  make  the  Eulerian  sectional  approach  attractive:  first,  the  sprays  used  for  mitigation  are 
generally  diffuse;  that  is,  the  number  density  and  size  of  particles  are  small  enough  that  the  volume  fraction 
of  water  is  very  low.  Because  of  this,  collisions  between  particles  have  a  minimal  influence  on  the  dynamics 
of  the  droplets,  and  are  ignored  for  this  particular  model.  Second,  although  the  volume  fraction  is  fairly 
small,  the  total  volume  of  the  domain  and  the  spray  region  can  be  quite  large,  thus  the  total  number  of 
droplets  or  particles  can  be  very  large.  Although  there  has  been  several  write-ups  describing  the  Eulerian 
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sectional  approach,  the  method  is  described  in  detail  below  because  of  differences  between  this  approach 
and  previous  sectional  approaches. 


3.1  Gas-Phase  Governing  Equations 

The  governing  equations  for  the  gas-phase  are  the  inviscid,  compressible  form  of  the  conservation  equa¬ 
tions  for  species,  mass,  momentum,  and  energy.  They  are  written  as: 


+  V-pkU  = 


-b  V  •  pu  =  5" 


+  V  •  puu  =  -  VP  +  5” 


^  +  V  ■iE  +  P)u  =  (7) 

where  5"*,  and  5®”*^  are  the  transfer  source  terms  between  the  dispersed-phase  and  the  gas- 

phase,  and  are  explained  below.  Given  N  species,  there  are  JV  -I-  3  conservation  equations  for  the  gas-phase 
flow  in  one-dimension.  For  the  simulations  completed  in  this  report,  the  gas-phase  is  assumed  to  be  an  ideal 
and  calorically  perfect  gas  mixture,  meaning  that  each  species  has  a  constant  specific  heat  and  obeys  the 
ideal  gas  law.  This  differs  from  our  previous  report  where  the  gas  was  only  assumed  to  be  thermally  perfect 
and  the  specific  heats  were  specified  as  functions  of  temperature.  The  relation  between  pressure,  species 
densities,  and  temperature  is  given  by  the  ideal  gas  law  as: 


P  =  RuT'£pk/Wk  (8) 

fc=i 

where  Rn  is  the  universal  ideal  gas  constant  (in  molar  units).  Because  we  have  a  calorically  perfect  gas,  the 
total  internal  energy  is  given  by  the  relation: 


^  1 
■^  =  X)  PkCv,kT  -h  -pu  •  u 


9 


The  specific  heat  for  each  species  is  specified  in  terms  of  the  specific  heat  ratio  (7*  =  Cp^k/Cy^k)  and  the 
ideal  gas  constant  {Rk  —  RulWk)  such  that  Cy^k  =  Rklil  —  !)•  We  use  a  multi-species  formulation  to 
differentiate  between  driver  gas,  driven  gas,  and  water  vapor.  For  air,  the  molecular  weight  is  Wair  =  28.97 
kg/kmol  and  the  specific  heat  ratio  is  jair  =  1-4,  and  for  water  vapor  the  molecular  weight  is  Wh^o  = 
18.015  kg/kmol  and  the  specific  heat  ratio  is  =  1-324. 

3.2  Dispersed-Phase  Governing  Equations 

In  a  very  general  form,  the  dispersed-phase  can  be  described  by  a  probability  density  function  at  the 
kinetic  level.  Assuming  the  spray  is  made  up  of  spherical  particles  characterized  by  a  single  geometry 
variable  ^  (which  generally  is  either  the  diameter,  surface  area,  or  volume  of  the  droplet),  single  velocity  u;, 
and  a  single  temperature  T/,  a  probability  density  function  /^(f,  x,  0,  u;,  Tj)  can  be  determined  such  that  the 
number  density  at  time  t  is  given  as  dt  dx d<l> dui  dTi  around  the  point  (t,  x,  /^,  rxuTi).  The  Boltzmann 
kinetic  transport  equation  for  the  probability  density  function  is  given  by: 

-f-  Vx  -  {uif^)  -t-  +  V„,  •  (F/^)  -I-  =  Tc  -f  Te  (10) 

where  ii^  is  the  rate  of  change  of  the  particle  geometry  (JR^  =  d<j)/dt)  due  to  vaporization,  F  is  the  momen¬ 
tum  exchange  between  the  dispersed-phase  and  gas-phase,  Q  is  the  energy  exchange  between  the  dispersed- 
phase  and  gas-phase  excluding  phase  change,  Fc  represents  coalescence  and  break-up  of  the  particles,  and 
Te  represents  elastic  collisions  between  the  particles.  For  this  research,  Tc  and  Fg  are  negligible  and  are  not 
considered  further.  All  of  these  functions  require  modeling  or  experiments  for  closure. 

From  this  general  level,  one  can  further  develop  the  equations  in  several  different  ways  to  generate  usable 
models  for  simulations.  The  method  we  choose  is  the  so-called  sectional  approach  originally  introduced  by 
Tambour  [9].  The  sectional  approach  divides  up  the  dispersed-phase  in  each  cell  into  M  discrete  sections 
based  on  the  size  of  the  particles.  Every  section  j  is  bounded  by  particle  sizes  and  Our  discussion 

of  the  spray  model  closely  follows  the  work  of  Laurent  and  Massot  [12].  Within  each  section,  we  make  the 
following  assumptions: 

1.  The  probability  density  function  of  the  spray  is  such  that  at  a  given  size  and  location  (x,  t,  <f>),  there  is 
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one  and  only  one  characteristic  velocity  and  temperature. 

2.  The  dispersion  of  the  distribution  function  at  a  given  size  is  zero  for  velocity  and  temperature. 

3.  In  one  section,  the  characteristic  velocity  and  temperature  do  not  depend  on  the  size  of  the  particles. 

4.  There  is  a  one-to-one  correspondence  between  surface  temperature  and  particle  enthalpy. 

For  each  section,  the  characteristic  velocity  and  temperature  are  denoted  and  T^.  These  assumptions  al¬ 
low  us  to  express  the  probability  density  function  in  a  much  simpler  form,  as  n(<,  x,  ^),  such  that  n(f,  x,  (^)  dt  dx  d<j) 
is  the  probable  number  of  particles  around  the  point  (<,x,  ^).  We  assume  that  the  liquid  or  solid  material 
is  incompressible,  thus  the  constant  pressure  and  constant  volume  specific  heats  are  equal.  The  specific 
internal  energy  of  the  particle  is  defined  as: 


Gd  —  CiTd  -h  Bf  i 

is  set  such  that  the  latent  heat  of  vaporization  Ly  is  defined  as: 


(11) 


Ly{T)  -  eH^0{g){T)  -  ejj20{i){T')  =  {Cy^H^oig)  -  Ci)T  -  (12) 

This  indicates  that  the  latent  heat  of  vaporization  should  vary  with  temperature,  depending  on  the  difference 
in  constant  specific  heats  of  the  vapor  and  liquid  phases. 

The  basic  mass,  momentum,  and  energy  equations  can  be  recovered  by  multiplying  Eqn.  10  by  1,  u^, 
and  Bd  and  integrating  over  velocity  and  temperature: 

dn  _  ,  .  d  ,  ^  . 

+  (nud)  +  =  0  (13) 

d  Q 

^(nud)  +  Vx  •  (nUrfUd)  -t-  —{nR^Ud)  =  nF  (14) 


d  .  s  „ 

— (ned)  -I-  Vj; 


(raUdCd)  -I-  —  (ni?d>6d) 


nCiQ 


(15) 
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For  the  above  equations  we  have  omitted  bars  associated  with  averaging  over  velocity  and  temperature 
because  of  the  above  assumptions  and  for  clarity.  To  obtain  the  conservation  equations  commonly  used  for 
the  sectional  approach,  we  integrate  with  respect  to  (f)  in  each  section.  The  number  density  of  the  section  is 
simply: 

(16) 

J(pU) 

Similarly,  one  can  get  the  mass,  momentum,  and  enthalpy  densities  of  the  dispersed-phase  by  integrating: 

r^O'+i) 

=  piV{(j))n{t,x,(f>)d^  (17) 

p^v  {(f>)n{t,  X,  ^)ud(f,  x)  #  (18) 

p^y  x,  4^)edit,  x)  #  (19) 

where  pi  is  the  density  and  V’(^)  is  the  volume  of  the  particles.  Other  terms  in  the  overall  conservation 
equations  are  similarly  treated: 


p^y X,  (j))F{t,  X,  (f>)  d<l)  (20) 

m^)QU)  =  piV{(l))n{t,x,  (f>)CiQ{t,  x,  <f))  d<f>  (21) 

...  ^  .N  Q 

mO)R(^)  =  p^yf^^)^  [n{t,x,(f>)Rd,{t,x,4>)]  d(f>  (22) 

Integration  by  parts  is  used  for  Eqn.  22  to  obtain  a  more  convenient  representation  of  the  vaporization: 

4(fl  '■‘''a#  -4„, 

This  allows  us  to  define  two  vaporization  terms,  given  below: 

=  -pin{t,x,<l>^^'>)R^{4>^^'>)V{(f)^^^)  (24) 


Physically,  represents  transport  of  mass  from  section  j  to  section  (j  -  1),  and  represents  transport 
of  mass  from  section  j  to  the  gas-phase. 
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The  conservation  equations  are  then  written  for  each  section  j  in  the  typical  sectional  form: 


dt 


(26) 


dt 


+  Va;  •  =  -  (e[^^  +  E^^'j  (28) 

+m(j)  (q(J)  -  E^Ly)  +  •  vif 


The  source  terms  for  the  gas-phase  are  simply: 


M 


i=i 


gmom  ^  ^  {EPm^^nP  -  F(^)m(j)) 

M 

genth  (EPm'^^ef  -  (q(^)  -  E^L^)  -  ■  xxf) 


(29) 

(30) 

(31) 


j=l 


The  species  source  term  is  simply  the  mass  source  term  for  the  water  vapor  species,  and  zero  for  everything 
else. 


3.3  Particle/Droplet  Models 

To  develop  expressions  for  E^,  E^,  F(^\  and  Q(j)  we  need  to  use  particle  and  droplet  models  to 
represent  the  processes  of  vaporization,  heat  transfer,  and  drag.  The  present  simulations  use  models  that 
have  been  used  widely  in  the  droplet  community.  We  describe  these  models  briefly  below.  The  models  are 
•s  based  on  drag  and  vaporization  around  a  single  particle  in  quiescent  flow,  that  are  then  corrected  for  different 

flow  regimes  based  on  empirical  correlations. 

Particle  drag  is  computed  based  on  Stokes  flow  with  high  velocity  correction  terms  [18].  For  a  single 
particle, 

=  37r/xur£>(l  -h  0.15ile°  ®®’')  (32) 

where  //  is  the  gas  viscosity  at  the  particle  temperature.  We  use  this  expression  to  obtain  an  expression  for 
F  found  in  Eqn.  10: 

F  =  ^Ur£>(l  -I-  0.15i?e°-®®'*’)  (33) 
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The  evolution  of  the  temperature  of  the  particle  is  given  as 

=  ^D^h[T  -  Trf]  +  pi^Ly  (34) 

where  h  =  NuX/D  is  the  heat  transfer  coefficient,  T  is  the  far-field  gas  temperature,  Ly  is  the  latent  heat 
of  vaporization,  is  the  particle  temperature,  and  Nu  is  the  Nusselt  number.  The  Q  in  Eqn.  10  is  the  heat 
transfer  from  the  gas  phase  into  the  liquid  phase,  and  does  not  account  for  phase  change,  thus 

The  Nusselt  number  can  either  be  assumed  to  be  2  (appropriate  for  non-convective  environments),  or  for 
convective  environments  it  can  be  represented  by  the  formula: 

Nu  =  2  +  0.459Pr^/®i2e°'®^  (36) 

The  vaporization  models  used  for  all  of  the  water  droplet  simulations  in  this  report  are  based  on  spherically- 
symmetric  isolated  droplet  solutions,  with  modifications  added  based  on  experimental  correlations  to  ac¬ 
count  for  convective  effects.  The  first  model  we  use  is  the  basic  Z)^-law  given  by: 

(37) 

where  D  is  the  particle  diameter  and  ^  is  the  vaporization  coefficient.  The  model  used  assumes  a  simple 
temperature  dependent  expression  for  the  vaporization  such  that: 

p{T)  =  /0o  [l  +  7.4233  x  10"'^  (T  -  300)^’^®'‘®]  (38) 

where  /0o  =  7600  /xm^/s,  obtained  from  [19,  20].  This  will  be  referred  to  as  the  empirical-/?  model,  because 
it  is  based  totally  on  empirical  data.  Although  a  simple  and  efficient  model,  a  significant  problem  with 
this  model  is  that  it  does  not  account  for  droplet  heat-up.  For  our  simulations,  as  the  shock  wave  passes  a 
given  droplet,  initially  the  droplet  vaporization  will  be  small  as  most  of  the  available  heat  is  used  to  heat 
up  the  droplet  to  near  its  boiling  point.  Only  after  the  droplet  is  heated  somewhat  will  vaporization  become 
significant.  This  heat-up  period  is  critical,  since  to  be  effective  the  droplet  must  vaporize  close  to  the  shock 
wave. 
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The  second  vaporization  model  addresses  this  issue.  It  is  based  on  the  classical  spherically-symmetric 
droplet  vaporization  model  of  Spalding  [21,  22],  but  includes  effects  from  droplet  heating.  The  model 
assumes  the  gaseous-H20  pressure  at  the  surface  of  the  droplet  is  the  equilibrium  vapor  pressure  of  water  at 
the  local  temperature,  and  is  thus  called  the  equilibrium  model  throughout  this  report.  Note  that  this  is  yet 
another  equilibrium  that  only  pertains  to  the  amount  of  water  vapor  at  the  droplet  surface.  It  is  not  in  any 
way  associated  with  the  gas-particle  equilibrium  or  the  shock  structure  equilibrium  mentioned  previously. 
Assuming  a  unity  Lewis  number,  the  vaporization  coefficient  P  is  calculated  as: 


O  \ 

P=—HB  +  l)  (39) 

Pl^P 

where  A  is  the  thermal  conductivity  of  the  gas,  pi  is  the  liquid  density,  Cp  is  the  constant  pressure  specific 
heat  of  the  gas,  and  B  is  the  Spalding  transfer  number,  defined  as: 

B  =  {YH20,g  -  Yh20,s)/ {Yh20,s  “  1)  (40) 

where  YH^o,g  is  the  mass  fraction  of  water  in  the  far-field,  Yff^o,s  is  the  mass  fraction  of  water  at  the  surface 
of  the  droplet.  Yji^o^s  is  calculated  from  the  partial  pressure  of  water  at  the  droplet  surface, 

_  _  PH20,sat{Td)WH20 

^H20,s  - - ^ -  (41) 

where  PH20,sat  is  the  saturation  pressure  of  water,  Ta  is  the  droplet  temperature,  WH2O  is  the  molecular 
weight  of  water,  and  W  is  the  molecular  weight  of  the  gas  mixture. 

Assuming  the  I?^-law,  the  spherically-synunetric  mass  vaporization  rate  for  a  single  droplet  is  calculated 


dV  ttD  dD'^  ttDB 

=  =  «-—  =  -«—  (42) 

where  rhss  represents  the  spherically  symmetric  solution.  For  convective  environments,  the  correlation  of 
Ranz  and  Marshall  [23]  can  be  used  to  determine  the  correct  vaporization  rate.  This  is  expressed  as: 

m  =  rhss  [l  +  O.ZPr^/^{2Re^/^)]  (43) 

where  V  is  the  droplet  volume,  Pr  =  pCp/X  is  the  gas-phase  Prandtl  number  at  the  droplet  temperature, 
and  Re  =  purD/ p  is  the  Reynolds  niunber  based  on  droplet  diameter,  gas  density  and  viscosity,  and  relative 
velocity  Ur  —  |ur|  =  |u  —  Ud|.  Note  that  all  unsubscripted  variables  are  gas  quantities. 
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For  the  vaporization  models,  the  heat  transfer  is  modified  by  an  additional  term  that  is  typically  expressed 


in  the  Nusselt  Number: 


Nu  =  Nu* 


ln(l  +  B) 
B 


(44) 


where  Nu*  is  the  original  Nusselt  number  for  non-vaporizing  particles  expressed  in  Eqn.  36.  The  additional 
term  does  not  change  the  expressions  developed  for  the  sectional  source  terms,  as  the  additional  term  is 
independent  of  droplet  geometry. 

For  the  drag  and  heat  transfer  computations,  the  Sutherland  expression  is  used  for  expressing  viscosity 
as  a  function  of  temperature  in  air: 


•V 


=  (45) 

where  /  =  1 .458  x  10~®  kg/(m  sec)  for  a  mixture  containing  predominantly  air.  The  thermal  conductivity 

is  calculated  assuming  a  constant  Prandtl  number  of  0.75, 


A(r)  = 


Pr 


(46) 


For  glass  particles,  the  only  physical  properties  we  need  are  the  density  of  glass  Pgiaas  =  2500  kg/m® 
and  the  constant  specific  heat  for  glass  Cgiass  =  766  J/kg  K.  For  liquid  water  droplets  has  a  constant 
density  of  Ph20{1)  =  1000  kg/m®  and  a  constant  specific  heat  of  Ch^o(1)  =  4184  J/kg  K.  The  latent  heat  of 
vaporization  is  =  2.3049  x  10®  J/kg  at  300  K.  The  saturation  pressure  of  water  PH20,sat{T)  is  curve  fit 
with  widely  available  data  and  has  a  form  similar  to  the  Clausius-Clapeyron  equation: 


PH20,aatiT)  =  A  exp  (47) 

where  for  water,  A  =  1.20967  x  10®  bar,  B  =  3835.83  K,  and  C  =  -45  K.  As  with  the  glass  particles, 
we  use  the  Sutherland  expressions  to  determine  the  viscosity  and  thermal  conductivity  as  a  functions  of  the 
temperature  only  (Eqns.  45  and  46). 

For  the  simulations,  a  cutoff  temperature  is  arbitrarily  chosen  below  which  the  vaporization  rate  becomes 
zero.  This  temperature  is  selected  above  3(X)  K  such  that  no  vaporization  takes  place  before  the  shock  wave 
reaches  the  droplets.  For  these  simulations,  a  value  of  Tc  =  310  K  was  chosen.  Lxiwer  values  between  3(X) 
and  310  K  could  also  be  chosen,  but  only  have  a  minor  effect  on  the  solutions. 
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3.4  Cross-coupling  Source  Terms  for  the  Eulerian  Sectional  Approach 

Given  models  that  describe  single-particle  behavior,  we  develop  expressions  for  the  source  terms  found 
in  the  sectional  conservation  Eqns.  26-29.  For  the  model  development,  we  choose  the  geometry  variable  (j) 
to  be  the  volume  V.  An  important  point  to  mention  is  that  the  resulting  source  terms  are  different  depending 
on  whether  (f)  represents  the  diameter  D,  surface  area  S,  or  volume  V.  In  the  work  of  Greenberg  and 
Tambour[23],  they  choose  V  as  their  geometric  variable,  however,  Laurent  and  Massot[12]  have  chosen  the 
surface  area  S  as  their  geometric  variable. 

The  fundamental  approximation  of  the  sectional  approach  assumes  a  relationship  between  the  number 
probability  density  function  n{t,  x,  (^)  and  the  geometry  variable  (j)  within  each  section.  The  most  reasonable 
choice  is  to  take  0)  to  be  a  constant  in  each  section.  This  simplifies  the  resulting  integration 
and  gives  a  clear  method  of  refinement.  Either  the  number  density  or  the  mass  density  formulas  can  be  used 
to  fix  the  value  of  n{t,  x,  for  section  j.  We  find  that  fixing  it  with  respect  to  the  mass  density  is 

the  most  convenient. 

Both  the  surface  area  S  and  the  particle  diameter  D  can  be  related  to  the  particle  volume  V  through: 

(48) 

(49) 

This  is  done  by  using  the  definition  of 

^yW+i) 

=  /  piVn{t,-K,V)  dV 

Jv(3) 

Noting  that  n  =  and  is  a  constant  with  respect  to  V,  this  equation  is  integrated  to  get 

solving  for  we  get 

Jj)  ^  2m^)/pi 

(y(i+i))2  _  (v0))2 


(50) 


(51) 


(52) 


For  calculating  the  vaporization  rate  terms  and  we  develop  an  expression  for  =  dV/dt 
based  on  the  £>^-law  (Eqn.  37).  Knowing  that  V  =  {it/6)D^  and  through  differentiation  and  algebraic 
manipulation  of  the  r>^-law,  we  obtain  an  expression  for  the  change  in  volume  for  the  spherically  synunetric 
droplet: 

^  (53) 

Using  the  droplet  vaporization  correlation  of  Ranz  and  Marshall  (Eqn.  43),  we  obtain  the  change  in  volume 
for  the  droplet  in  a  convective  environment  as: 

^  ^  [l  +  0.3Pri/3(2i2e^/2)j  ^  ^  ^  0.3Pr^/3(2ile^/2)J  (54) 


Pulling  the  droplet  diameter  out  of  the  Reynolds  number,  substituting  this  and  the  value  for  into 
Eqn.  24  we  get  an  expression  for  e[^^  as 

where  A  =  Pr^/^{2pUr / .  By  a  similar  means,  we  substitute  into  the  definition  for  E^^  to 


obtain  the  integral: 


El 


(i) 


yO+i) 


rV'~‘ 


PlK 


W 


[1  +  0.3Pr^/^(2i?e^/2)]  dV 


(56) 


Integrating  this  and  evaluating  the  result  at  and  V^\  we  obtain: 


"T 


(y(j+i))2  _  (u(i))2 


[(y(j+l))4/3  _  (^^(i))4/3j  4.  2^  [(y(j+l))3/2  _  (y(i))3/2j  1 


J 

(57) 


Using  the  particle  drag  model  (Eqn.  32),  we  develop  an  equation  for  the  integrated  force  felt  on  section 
j  as  defined  in  Eqn.  20.  The  function  F  found  in  Eqn.  20  is  the  acceleration  felt  by  a  particle  due  to  drag, 
and  is  thus: 


F(F) 


37r//U; 


r 


PiV 


1+0. 


Uj  {-V 


0.229 


(58) 


To  solve  for  the  cumulative  force,  we  take  Eqn.  20  and  substitute  the  above  expression  in  for  F,  substitute 
«;(•’)  for  n,  and  substitute  U  for  P  =  ^GU/tt  to  get 


=  C"  0  +  (1)“"'  (t 


dV 


(59) 
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This  temi  is  then  split  into  two  terms,  and  for  simplicity.  Integration  of  the  F^  equation  with 
respect  to  V  is  easily  done,  and  with  some  rearrangement  and  substitution  of  we  get  the  following 


expression: 


,fa)FM  =  r(yO-«))4/3-(y(J))4/s1 


^  Pi)!  4  [  (yo+i))2  -  (y(j))2 

(i) 

Detenmnation  of  the  F^^^  is  also  simple,  although  the  algebra  is  a  bit  messier.  First,  we  substitute  in  for  the 
Reynolds  number  Re/D  =  pUr/p,  and  integrate  with  respect  to  volume: 

^  [{ye+»)i«a  - (fW))1.442s]  (61) 

substituting  in  for  k^^'>  and  multiplying  out  the  initial  numbers  we  get  the  expression: 


(j)  _  (,-^2.604  r (F(J+1))1-5623  _  (y(j))1.5623 


(f)' 


PI  [  (F(i+i))2-(F(i))2  jv^y 

The  last  expression  we  need  is  for  For  this,  we  need  an  expression  for  Q,  the  change  in 

temperature  with  respect  to  time.  An  equation  is  developed  from  Eqn.  34  as: 


Expressing  this  in  terms  of  the  volume  of  the  particle  V, 


0.55  /ai/\  0.183 


(f) 

ilco  CT\lit  tViic  iTitrx  txx/rk  tpninc  nO)  onrl  nU) 


This  expression  is  also  split  this  into  two  terms,  QY  and  to  simplify  the  integrals: 


Jv(i)  2Ci 


Doing  the  integration,  and  substituting  in  for 


,WqW  ^  r(l/W-))4/3  _  (l/(,))4/3l 

4Cipi  (fo+i))2  -  (i/a))2 


The  second  term  can  be  obtained  similarly, 


yO)  Cl  \7CJ  \  P  J 
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Again,  by  integrating  and  substituting  into  we  obtain  a  more  manageable  expression: 

PiCi  \  n  J 

The  above  expressions  for  and  close  the  system  of  conservation  Eqns.  26  -29,  and 

allows  us  to  solve  the  system  of  equations. 

4  Solution  Procedure 

The  solution  procedure  involves  solving  M  +  1  sets  of  conservation  equations,  where  M  is  the  number 
of  sections  being  simulated.  Each  set  contains  a  continuity  equation,  momentum  equations,  and  an  energy 
equation.  The  sets  are  coupled  together  through  the  source  terms  described  above.  The  present  procedure 
used  for  solving  these  conservation  equations  is  by  a  time-step  splitting  method.  Each  set  of  conservation 
equations  is  solved  independently  without  the  cross-coupling  source  terms.  Because  the  gas-phase  and 
dispersed-phase  can  have  significantly  different  time-step  restrictions,  the  gas-phase  solution  is  allowed  to 
take  several  time-steps  for  each  dispersed-phase  time-step  taken.  The  cross-coupling  source  terms  are  added 
at  the  end  of  each  dispersed-phase  time-step. 

4.1  Gas-phase  solution  procedure 

The  solution  procedure  used  for  the  gas-phase  conservation  Eqns.  5-7  is  the  FCT-algorithm  of  Boris 
and  Book  [24],  which  is  especially  well  suited  for  high-speed  flows.  The  version  of  the  algorithm  used  in 
this  project  is  described  in  detail  in  NRL/MR/64 10-93-7 192  [25],  and  will  not  be  repeated  here  for  brevity. 
Unlike  the  computations  described  in  the  memorandum  report,  we  use  the  explicit  ideal  gas  law  (Eqn.  8) 
and  the  energy  definition  (Eqn.  9)  instead  of  combining  the  two  to  obtain  an  expression  for  pressure  based 
on  internal  energy. 

4.2  Dispersed-phase  solution  procedure 

The  solution  procedure  used  for  the  dispersed-phase  conservation  Eqns.  26-29,  is  also  based  on  the 
explicit-FCT  algorithm  used  above.  Examining  the  conservation  equations,  notice  that  the  momentum  equa¬ 
tion  is  decoupled  from  the  energy  equation.  One  benefit  of  this  is  that  the  Courant  number  restriction  for  the 


1.5167 


(yW) 


1.5167 


(y(i+l))2  _  (^t3))2 


XTr 


(68) 
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numerical  procedure  is  based  on  the  particle  velocity  instead  of  the  acoustic  velocities  as  for  the  gas- 
phase,  allowing  for  much  larger  time-steps.  The  internal  energy  and  the  temperature  of  the  dispersed-phase 
is  similarly  related  through  a  constant  specific  heat  +  e°f  {). 

4.3  Gas/dispersed-phase  coupling 

The  gas/dispersed-phase  coupling  terms  are  only  added  to  the  solution  after  each  particle-section  ad¬ 
vances  one  time-step,  and  the  gas-phase  has  advanced  multiple  time-steps.  The  coupling  terms  are  then 
added  explicitly.  Each  section  of  the  dispersed-phase  is  coupled  to  the  neighboring  sections  through  the 
vaponzation  term  EY'  (Eqn.  24).  Every  section  is  coupled  to  the  gas-phase  through  the  second  vaporization 
term  (Eqn.  25),  the  momentum  transfer  term  (Eqn.  20),  and  the  heat  transfer  term  (Eqn.  21). 
For  our  selection  of  <p  =  V  and  with  the  droplet/particle  models  presented  above,  the  final  expression  for 
Ei^  is  given  in  Eqn.  55,  E^^  is  given  in  Eqn.  57,  the  drag  somce  term  is  given  in  Eqns.  60  and  62,  and 
the  heat  transfer  source  term  is  given  in  Eqns.  66  and  68.  The  complete  source  terms  introduced  to  the 
gas-phase  conservation  equations  are  given  in  Eqns.  29-31. 

4.4  Initial  and  Boundary  Conditions 

The  simulations  for  this  report  are  all  one-dimensional  shock  tube  simulations.  The  basic  geometry 
is  shown  in  Figure  3.  The  solution  domain  is  divided  into  three  sections:  a  driver  section,  an  unseeded 
driven  section,  and  a  seeded  driven  section.  In  each  section  the  temperature,  pressure,  and  composition  are 
uniform  initially.  The  seeded  region  contains  uniformly  distributed  particles  at  a  specified  mass-loading. 
The  velocity  for  the  entire  domain  is  initially  zero,  and  the  particles  are  assumed  to  be  perfectly  suspended 
in  the  gas  with  zero  velocity.  This  type  of  simulation  is  equivalent  to  experiments  involving  shock  tubes 
with  diaphragms,  where  a  diaphragm  initially  separates  the  driver  and  driven  gases  and  is  burst  at  the  start 
of  the  experiment.  The  driver  gas  is  at  a  much  higher  pressure,  and  thus  pushes  the  driven  gases  much  like 
a  piston,  creating  a  shock  that  propagates  through  the  driven  gases. 

Boundary  conditions  for  each  end  of  the  shock  tube  are  perfectly  reflecting  walls.  Typically,  the  sim¬ 
ulation  is  stopped  before  the  shock  wave  reflects  off  of  the  right  wall.  The  rarefaction  waves  do,  however, 
reflect  off  of  the  left  wall  and  can  subsequently  influence  the  shock  wave. 
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5  Numerical  Model  Validation 


Numerical  validation  of  the  computational  scheme  used  the  exact  Sommerfeld  geometry  [7],  consisting 
of  a  2  m  long  driver  section,  a  1.05  m  long  unseeded  driven  section,  and  a  4.76  m  long  seeded  driven  section. 
The  driven  section  was  seeded  with  mono-dispersed  27  /xm  glass  particles.  Several  experiments  were  run 
with  varying  amounts  of  mass-loading,  and  with  different  driver  pressures.  The  driver  section  pressures  are 
indicated  below  for  each  case.  The  resolution  of  the  simulation  is  approximately  2.5  mm  per  cell,  enough  to 
obtain  accurate  solutions,  but  it  does  cause  some  noise  in  the  shock  wave  Mach  number  calculation. 

Sommerfeld  was  interested  primarily  in  the  transition  region  between  a  shock  wave  traveling  through 
unseeded  flow  and  seeded  flow,  as  is  shown  in  Figure  4.  In  the  transition  region,  the  shock  wave  deceler¬ 
ates  from  the  unseeded  shock  wave  velocity  to  the  equilibrium  shock  wave  velocity  in  the  seeded  mixture. 
Because  the  energy  and  momentum  transfer  are  not  instantaneous  between  the  gas  and  dispersed-phase,  the 
shape  of  the  pressure  and  temperature  profiles  behind  the  shock  wave  also  changes  due  to  the  particles.  To 
ensure  the  expansion  fan  did  not  reflect  back  and  impact  the  front  shock,  he  used  a  fairly  long  driver  section 
(which  is  typical  of  many  shock  tube  experiments).  In  addition,  Sommerfeld  used  a  fairly  large  unseeded 
driven  section  to  ensure  that  the  shock  wave  was  properly  formed  before  impacting  the  seeded  gas. 

In  Figure  5,  we  show  comparisons  of  shock  wave  Mach  number  between  four  separate  Sommerfeld 
experiments  and  their  corresponding  simulations.  The  comparison  for  all  four  cases  is  fairly  good.  The 
experimental  results  (shown  as  symbols  on  each  of  the  four  plots),  shows  that  the  transition  region  can  be 
very  roughly  divided  into  three  almost  linear  segments.  Initially,  the  deceleration  is  very  fast.  At  about 
one  meter  past  the  seeded-unseeded  interface,  the  deceleration  abraptly  slows.  At  about  three  meters  past 
the  interface,  the  deceleration  again  slows  as  it  approaches  the  shock  structure  equilibrium  value  for  the 
gas-particle  mixture  (found  from  the  equivalent-gas  analysis  presented  earlier). 

For  the  simulations,  the  initial  deceleration  is  accurately  captured  for  all  cases.  The  abrupt  change  in  the 
deceleration  is  not  totally  caught  in  all  the  numerical  simulations,  resulting  in  slightly  lower  Mach  numbers 
in  the  calculations  as  compared  to  the  experiments  for  the  Mach  number  1.49  and  1.26  cases  with  mass¬ 
loading  of  0.63.  Interestingly,  the  higher  mass-loading  cases  and  higher  velocity  cases  capture  this  shift  in 
deceleration  fairly  accurately,  and  do  not  show  the  same  type  of  errors.  The  poorest  comparison  is  with  the 
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Figure  5:  Comparison  between  numerical  simulations  and  experiments  [7]  for  shock  wave  Mach  numbers. 
The  air  is  seeded  with  27  ^m  glass  particles  with  the  given  mass-loading. 


low  Mach  number  case  (Fig.  5d),  although  there  is  no  clear  reason  why  this  discrepancy  should  exist  for  this 
case.  Our  results  are  consistent  with  other  numerical  simulations  of  the  Sommerfeld  simulations,  presented 
by  Sommerfeld  and  other  researchers,  indicatiug  that  the  problem  is  either  with  the  underlying  model  for 
low  Mach  number  shock  waves  or  with  the  experiment  itself,  and  not  with  our  implementation  of  the  model. 

The  effect  of  including  particle  heating  in  the  model  is  shown  clearly  in  Figure  6.  Within  two  meters 
of  the  start  of  the  seeded  section,  the  simulation  results  are  very  close,  with  the  simulation  without  particle 
heating  actually  comparing  slightly  better  with  experiments  than  the  simulations  with  particle  heating.  Fur¬ 
ther  downstream,  however,  the  case  without  particle  heating  underpredicts  the  attenuation.  This  is  because 


Figure  6:  Effect  of  particle  heating  on  shock  wave  Mach  number.  Mo  =  1.49, 77  =  0.63,  26-28  fim  glass 
particles. 

of  the  additional  energy  that  is  taken  out  of  the  gas  flow  to  equilibrate  the  temperatures  for  the  simulations 
with  particle  heating.  This  energy  is  never  removed  for  the  case  without  particle  heating,  resulting  in  an 
incorrect  “equilibrium”  shock  wave  Mach  number. 

A  possible  reason  for  the  better  comparison  early  in  the  transition  region  for  the  simulation  without 
heat  transfer  could  be  because  of  the  infinite  conductivity  model  used  for  the  particle  heating.  This  model 
keeps  the  surface  temperature  of  the  particles  artificially  low,  increasing  the  heat  transfer  and  decelerating 
the  shock  wave  prematurely.  As  the  particles  start  to  equilibrate,  the  infinite  conductivity  model  becomes  a 
better  approximation,  and  the  experimental  and  numerical  results  start  to  show  better  agreement.  One  of  the 
limitations  of  the  sectional  approach  is  the  necessity  of  using  an  infinite  conductivity  approach,  instead  of 
more  complex  temperature  models  for  the  droplets.  This  issue  also  appears  when  discussing  the  equilibrium 
vaporization  model  below. 

In  the  original  Sommerfeld  paper,  the  experimental  procedure  for  producing  a  pressure  ratio  of  four 
between  the  driver  and  driven  gas  is  not  specified.  In  all  the  simulations  to  this  point,  an  assumption  was 
made  that  the  pressure  was  raised  through  heating  the  driver  gas  (and  having  the  density  remain  constant). 
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Figure  7;  Comparison  of  unmitigated  and  mitigated  shock  wave  Mach  numbers.  The  initial  pressure  differ¬ 
ence  for  the  shock  tubes  is  obtained  through  temperature  =  4,  P4//91  =  1)  or  density  {paIp\  =  4, 

T4/Ti  =  1).  P4/P1  =  4, 77  =  0.63.  26-28  pm  glass  particles. 

Another  procedure  to  obtain  an  initial  pressure  ratio  of  four  is  to  compress  the  gas  and  keep  the  temperature 
constant.  These  two  different  initial  conditions  are  compared  below  in  Figures  7  and  8  with  (dashed)  and 
without  (solid)  particle  seeding.  There  is  considerable  difference  between  the  solutions  due  to  the  initial 
conditions,  although  the  transition  behavior  is  very  similar  between  the  two  cases.  From  the  results  above 
it  is  apparent  that  Sommerfeld  used  heating  to  obtain  the  correct  pressure  difference  between  the  driver  and 
driven  gases.  However,  we  present  both  results  for  completeness,  and  note  that  most  shock  tubes  generally 
use  a  combination  of  heating  and  compression  to  obtain  the  desired  pressure  difference  between  driver  and 
driven  gases. 

Some  insight  into  the  reasons  for  the  difference  are  apparent  from  examining  pressure  profiles  as  shown 
in  Figure  8.  The  left  side  corresponds  to  an  initial  driver  temperature  ratio  of  T4/T1  =  1  and  a  density 
ratio  of  pi/pi  =  4,  and  the  right  side  corresponds  to  the  inverse  =  4,  p^/p\  =  1).  Note  that  the 

second  case  is  simply  the  validation  case  presented  above  in  Figure  5b.  Although  both  simulations  used  the 
same  driver  pressure,  the  shock  pressure  that  develops  for  each  case  is  significantly  different,  resulting  in 
the  different  shock  wave  Mach  numbers  observed  in  Figure  7.  The  dashed  lines  show  computations  with 
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Figure  8:  Snapshot  of  solutions  for  unmitigated  (solid)  and  mitigated  (dashed)  shocks.  The  initial  pressure 
difference  for  the  shock  tube  is  obtained  through  temperature  (a)  or  density  (b).  P4/P1  =  4, 7;  =  0.63, 
26-28  nm.  glass  particles. 


particle-seeding  of  the  flow,  and  clearly  show  that  in  both  cases  the  shock  wave  slows  down  after  it  enters 
the  seeded  region. 

The  Sommerfeld  experiments  (and  many  other  experiments)  focused  on  the  shock  wave  Mach  number 
in  the  transition  region  of  the  flow-field.  Although  fliis  is  an  important  aspect  of  mitigation,  also  of  interest  is 
the  shock  overpressure  and  maximum  overpressure  felt  at  any  given  point  in  the  domain.  While  for  unseeded 
shocks  (and  shocks  seeded  with  infinitely  small  particles),  shock  overpressure  and  maximum  overpressure 
are  equivalent,  for  gases  seeded  with  finite-size  particles  the  maximum  overpressure  occurs  downstream  of 
the  shock  front  (shown  in  Figure  8  and  discussed  below).  In  particular,  we  are  interested  in  the  pressure 
response  away  from  the  diaphragm  in  the  decay  region  of  the  shock. 

Although  the  Sommerfeld  geometry  is  useful  for  comparison  with  experimental  results,  it  is  not  the  best 
geometry  to  use  for  gaining  a  complete  understanding  of  characteristics  of  shock-particle  interactions.  There 
are  several  reasons  for  this.  Although  the  driver  section  is  long  enough  such  that  the  expansion  fan  does  not 
effect  the  shock  firont  and  shock  wave  Mach  number  (the  quantity  that  Sommerfeld  was  most  interested 
in),  it  does  effect  the  area  just  behind  the  shock  wave,  making  it  more  difficult  to  interpret  temperature  and 
pressure  profiles  obtained  at  different  time  intervals.  This  is  especially  apparent  for  the  droplet  vaporization 
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cases.  In  addition,  the  length  of  the  driven  section  is  such  that  only  the  transition  region  for  the  smallest 
particles  and  droplets  can  be  captured.  Finally,  because  we  are  also  interested  in  the  decay  region  of  a  shock 
wave  produced  from  a  shock  tube,  some  simulations  require  a  much  smaller  driver  section  and  much  longer 
driven  section.  Because  of  this,  our  results  are  split  into  two  sections.  One  section  describes  results  with  a 
very  large  driver  section.  This  is  done  to  understand  different  aspects  of  shock  waves  impinging  on  seeded 
flow.  The  second  set  of  results  examines  the  decay  region  of  the  shock  wave  produced  from  a  shock  tube. 
All  of  the  simulations  for  these  results  have  a  short  driver  section  and  long  driven  sections,  and  examine 
how  the  pressure  decays  as  the  shock  wave  propagates  through  the  driven  section.  Because  of  our  interest  in 
mitigating  the  shock  wave  with  either  glass  particles  or  water  droplets,  the  seeded  simulations  are  referred 
to  as  the  mitigated  cases,  and  the  unseeded  simulations  as  the  unmitigated  cases. 

6  Propagation  of  Shock  Waves  through  Seeded  Gases 

The  geometry  used  for  these  simulations  has  a  very  large  driver  section  (3  m),  a  small  unseeded  driven 
section  (0.1  m),  and  a  seeded  section  of  3-4  m.  The  small  size  for  the  unseeded  driven  section  is  adequate 
for  the  one-dimensional  simulations  to  form  the  shock  wave  before  it  impinges  on  the  seeded  flow.  Also, 
the  small  seeded  region  allows  us  to  examine  the  effect  of  large  particles  or  droplets  penetrating  into  the 
high-temperature  driver  gas.  The  resolution  used  for  these  simulations  was  1  mm  per  cell,  which  is  good 
enough  to  give  a  smooth  shock  wave  Mach  number  variation.  Simulations  computed  with  this  geometry  are 
indicated  in  Table  1.  Only  selected  simulations  are  discussed  in  detail. 

6.1  Representative  case  with  glass  particles 

A  representative  solution  (corresponding  to  case  #1  in  Table  1)  is  first  examined  with  the  long  driver 
section  mentioned  above.  The  seeded  regions  contains  a  nearly  mono-dispersed  flow  with  one  section  and 
a  diameter  range  from  26  fjim.  to  28  /zm.  The  mass-loading  for  all  of  these  cases  is  jy  =  =  0.63. 

Snapshots  of  the  solution  for  this  specific  case  are  shown  in  Figure  9  and  Figure  10.  Solution  snapshots  were 
taken  at  1  ms  intervals,  to  show  the  development  of  the  shock  wave  as  it  interacts  with  the  glass  particles. 
Figure  9  shows  the  unmitigated  shock  wave  results  (solid  lines)  compared  with  the  mitigated  shock  wave 
results  (dashed  lines).  Figure  10  shows  the  gas-phase  results  (solid  lines)  compared  with  the  dispersed-phase 
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Table  1:  Simulations  computed  with  the  long  driver  section  geometry.  Length  of  driver  section  is  3  m,  length 
of  unseeded  driven  section  is  0.1  m,  and  length  of  driven  section  is  3  m. _ 


Case# 

Dispersed-Phase 

Mass 

Loading 

Number 

Sections 

Particle 
Size  (/im) 

Vaporization 

Model 

1 

None 

0 

NA 

NA 

2 

Glass 

0.63 

1 

26-28 

2 

Glass 

0.63 

1 

13-14 

3 

Glass 

0.63 

1 

6.5-7 

4 

Glass 

0.63 

NA 

0 

5 

Glass 

0.63 

1 

4-6 

6 

Glass 

0.63 

1 

14-16 

7 

Glass 

0.63 

1 

24-26 

8 

Glass 

0.63 

1 

34-36 

9 

Glass 

0.63 

1 

49-51 

10 

Glass 

0.63 

1 

99-101 

11 

Water 

0.63 

6 

0-30 

None 

12 

Water 

0.63 

6 

0-30 

Empirical-;0 

13 

Water 

0.63 

6 

0-30 

Equilibrium 

14 

Water 

0.63 

6 

0-60 

Empirical-)8 

15 

Water 

0.63 

6 

0-15 

Empirical-y3 

16 

Water 

0.63 

6 

0-120 

Empirical-)9 

17 

Water 

0.63 

12 

0-30 

Empirical-;0 

18 

Water 

1.26 

6 

0-51.3 

Empirical-;0 

19 

Water 

0.315 

6 

0-17.54 

Empirical-;0 

20 

Water 

0.63 

6 

0-60 

Equilibrium 

21 

Water 

0.63 

6 

0-15 

Equilibrium 

22 

Water 

0.63 

6 

0-120 

Equilibrium 

23 

Water 

0.63 

12 

0-30 

Equilibrium 

24 

Water 

1.26 

6 

0-51.3 

Equilibrium 

25 

Water 

0.315 

6 

0-17.54 

Equilibrium 
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results  (dashed.lines)  for  the  mitigated  case. 

The  glass  particles  tend  to  have  the  maximum  effect  on  the  gas  flow  at  the  very  front  of  the  shock 
wave,  where  the  velocity  and  temperature  differences  between  the  gas  and  the  beads  are  largest.  Behind  the 
shock  wave,  the  velocity  and  temperature  of  the  beads  (and  gas)  adjusts  such  that  the  momentum  and  energy 
exchange  is  lessened.  Because  of  this  effect,  the  solution  snapshots  appear  rounded  at  the  front  of  the  shock 
wave,  instead  of  having  flat  profiles  as  they  do  for  the  unmitigated  cases  (this  is  clearly  shown  in  Figure  9.) 
Note  that  there  still  is  a  shock  front  with  an  associated  pressure,  density,  and  temperature  jump,  but  the  jump 
is  much  smaller  than  for  the  unmitigated  case.  It  is  also  much  lower  than  the  maximum  pressure  developed 
behind  the  shock  waves  for  the  mitigated  cases. 

An  interesting  aspect  of  the  mitigation  studies  is  that  the  maximum  pressure  behind  the  shock  wave  is 
actually  greater  for  the  “mitigated”  cases  than  for  the  non-mitigating  cases  (Figure  9c).  The  main  reason 
for  this  seems  to  be  an  increase  in  gas  density  behind  the  shock  wave  for  the  mitigating  cases  (Figure 
9a).  This  effect  is  predicted  from  the  equivalent-gas  analysis  above,  and  is  due  to  the  compression  wave 
generated  by  the  gas  slowing  down.  The  left-moving  reflected  shock  wave  mentioned  in  the  equivalent-gas 
analysis  is  seen  in  the  1  ms  profiles  at  an  axial  location  of  aiound  0.4  m.  It  can  be  seen  most  easily  in  the 
velocity  profiles  in  Figure  9f.  Although  a  compressive  shock,  it  is  not  sharp  because  momentum  and  energy 
equilibration  is  not  instantaneous  between  the  particles  and  gas. 

The  dispersed-phase  solution  is  shown  compared  with  the  gas-phase  solution  in  Figiure  10.  Here  we  see 
that  the  dispersed-phase  eventually  penetrates  into  the  driver  gas,  as  shown  clearly  in  the  density  (Fig.  10a) 
and  temperature  (Fig.  10b)  plots.  The  mass-loading  of  the  dispersed-phase  tends  to  spike  to  a  maximuTn 
within  the  driver  gas.  This  is  caused  by  the  larger  viscosity  within  the  driver  gas  (due  to  the  higher  tem¬ 
peratures)  increasing  the  amount  of  drag  felt  by  the  glass  particles.  Whether  the  particles  penetrate  into  the 
driver  gas  depends  on  the  size  of  the  unseeded  section  and  the  size  of  the  particles,  as  we  show  later. 

An  X  —  t  diagram  of  the  location  of  the  shock  wave  (green),  the  contact  discontinuity  (red),  and  the 
particle  interface  (blue),  as  well  as  the  shock  wave  Mach  number  is  shown  in  Figure  11.  The  mitigated 
case  is  shown  as  a  dashed  line.  Here  we  can  clearly  see  that  the  glass  particles  penetrate  past  the  contact 
discontinuity,  and  that  both  the  contact  discontinuity  and  the  shock  wave  are  slowed  considerably  by  the 
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Figure  9:  Snapshot  of  representative  solution  at  1  ms  intervals  for  unmitigated  (solid)  and  mitigated  (dashed) 
shocks.  Mo  =  1.49,  rj  =  0.63, 26-28  fim  glass  particles  with  large  driver  section  geometry. 
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Figure  10:  Snapshot  of  representative  solution  at  1  ms  intervals  fo 
(dashed)  for  the  mitigating  case.  Mo  =  1.49,  r}  =  0.63,  26-28  /rm 
geometry. 
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Figure  —  t  diagram  (a)  and  shock  wave  Mach  number  (b)  for  representative  mitigated  and  unmitigated 

cases.  For  the  a:  —  f  diagram,  green  indicates  shock  wave  location,  red  indicates  contact  discontinuity 
location,  and  blue  indicates  particle  interface.  Mo  =  1.49,  r)  =  0.63, 26-28  /xm  particles  with  the  extended 
driver  section  geometry. 


particles.  For  the  representative  case,  the  shock  wave  Mach  number  does  not  reach  the  equilibrium  value 
within  the  domain  (Figure  1  lb). 

6.2  Glass  particle  size  effects 


In  addition  to  the  above  representative  case,  several  other  simulations  were  run  with  a  similar  geometry 
to  determine  the  effect  of  particle  size  and  mass-loading  on  the  shock  wave  characteristics.  Particle  size 
has  a  significant  impact  on  the  mitigation  produced  by  the  glass  particles.  Changing  the  size  of  the  glass 
particles  changes  the  relaxation  time  required  for  the  particles  to  reach  equilibrium  with  the  gas  in  terms  of 
temperature  and  velocity.  This  affects  how  well  the  particles  slow  the  initial  shock  wave  down  and  penetrate 
into  the  driver  gases.  We  examined  droplet  diameters  ranging  from  100  /xm  to  infinitely  small,  with  a  mass¬ 
loading  ofr}  =  0.63.  For  the  infinitely  small  particles,  we  assume  the  glass  particles  are  in  thermal  and 
velocity  equilibrium  with  the  gas  at  all  times. 

Results  from  the  simulations  are  shown  in  Figure  12  and  13.  For  infinitely  small  particles,  the  shock 
structure  remains  a  square  wave,  although  the  shock  wave  itself  is  slowed  down  due  to  the  additional  mass 
being  carried  along  with  the  shock.  Although  the  shock  pressure  varies  considerably  with  particle  size,  the 
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maximum  pressure  within  the  driven  gases  remains  essentially  the  same  and  is  independent  of  the  particle 
size.  Comparing  the  shock  front  pressure  profile  (Fig.  12a)  at  1  ms,  2  ms,  3  ms,  and  4  ms,  the  profile  for  the 
simulation  with  27  /im  particles  changes  as  the  shock  wave  moves  downstream.  The  shock  front  pressure 
profile  for  the  simulation  with  infinitely  small  particles  does  not  change,  and  the  profile  for  the  simulation 
with  very  small  particles  (6.75  ^m)  only  has  small  changes  as  the  shock  moves  downstream.  Again,  this  is 
related  to  the  relaxation  time  associated  with  the  different  droplet  size.  At  a  certain  point  downstream,  each 
profile  reaches  an  equilibrium  value  after  which  the  profile  does  not  change  (assuming  the  driver  section 
is  long  enough  that  rarefaction  waves  do  not  effect  the  shock  front)  and  the  shock  wave  Mach  number  is 
constant. 

Figure  13  shows  the  s  —  t  diagram  and  the  shock  wave  Mach  number  as  a  function  of  location  for 
similar  cases.  Note  that  both  the  5  jjtm  and  the  infinitely  small  particle  cases  shock  Mach  number  have 
settled  to  the  same  equilibrium  Mach  number.  Also  note  that  the  contact  discontinuity  has  slowed  down 
by  almost  exactly  the  same  amount  for  all  the  mitigation  cases  regardless  of  the  particle  size.  Note  that  the 
infinitely  small  particles  very  quickly  reach  the  equilibrium  value,  whereas  the  finite  size  particles  require 
considerable  distance  to  reach  this  equilibrium  value.  This  equilibrium  value  is  only  weakly  dependent  on 
the  particle  size,  and  is  much  more  dependent  on  the  mass-loading  of  particles. 

Only  the  largest  particles  considered  actually  penetrated  into  the  driver  section  for  these  cases.  The 
amount  that  the  particles  penetrate  into  the  driver  gas  is  seen  both  in  the  ar  —  f  diagram  (Fig.  13a)  and  in  the 
driver-gas  density  plot  (Fig.  12b).  Only  for  the  25  ^m  case  have  the  particles  penetrated  into  the  driver  gas. 

6.3  Representative  case  with  water  droplets 

Realistic  mitigation  systems  would  use  water  droplets  and  not  solid  particles,  as  in  the  simulations 
described  previously.  Our  first  results  using  water  droplets  for  mitigation  show  snapshots  of  the  solution 
with  no  vaporization,  empirical-;^  vaporization  (Eqn.  38),  and  with  the  equilibrium  vaporization  (Eqn.  39). 
The  density  profiles  are  plotted  in  Figure  14,  and  the  temperature,  pressure,  and  velocity  profiles  are  plotted 
in  Figure  15.  As  shown  in  Figure  14,  the  total  density  and  the  density  of  the  air  both  rise  sharply  for  the 
vaporization  cases  compared  to  the  non-vaporizing  case.  This  is  because  vaporization  adds  mass  to  the 
gas-phase  and  also  compresses  the  gas-phase  by  slowing  the  shock  wave  down  further.  Looking  at  the 
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Figure  12:  Snapshot  of  pressure  and  driver  density  solutions  at  1  ms  intervals  for  shock  waves  mitigated  with 
different  size  particles.  Infinitely  small  (solid),  6.75  //m  (dashed),  27  ftm  (dash-dot)  particles  are  compared 
witfi  the  no  mitigation  solution  (dotted).  Mq  =  1.49,  rj  =  0.63,  with  the  laige  driver  section  geometry. 


Figure  13:  a:  —  f  diagram  (a)  and  shock  wave  Mach  number  (b)  for  several  different  size  glass  particles.  For 
the  af  -  <  diagram,  green  indicates  shock  wave  location,  red  indicates  contact  discontinuity  location,  and 
blue  indicates  particle  interface.  Mo  =  1.49,  rj  =  0.63,  with  the  larger  driver  section  geometry. 
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water  vapor  density,  we  see  that  the  empirical-/?  vaporization  model  evaporates  considerably  more  of  the 
droplets  than  the  equilibrium  vaporization  model,  and  most  likely  represents  an  upper  bound  on  the  amount 
of  vaporization  taking  place. 

Some  interesting  aspects  of  Figure  15  should  be  pointed  out.  For  both  vaporization  cases,  the  tempera¬ 
ture  is  reduced  significantly.  For  the  empirical-/^  vaporization,  the  temperature  in  the  shocked  driven  gases 
goes  down  to  the  cut-off  temperature  of  310  K,  at  which  point  no  further  vaporization  takes  place  and  the 
temperature  remains  constant.  Although  the  temperature  and  density  have  both  changed  significantly  with 
vaporization,  the  effects  of  these  changes  tend  to  cancel  out  such  that  the  resulting  shock  pressure  remains 
approximately  the  same  (it  is  actually  slightly  lower  for  high  levels  of  vaporization). 

Another  interesting  result  is  the  total  water  mass  in  each  of  the  sections.  This  quantity  tells  us  how  much 
water  has  evaporated,  and  how  much  mass  has  transferred  between  the  different  sections.  The  integrated 
mass  is  shown  in  Figure  16.  Because  of  the  large  differences  in  mass  between  sections,  the  y-axis  is  shown 
as  a  log  plot,  although  within  each  section  the  relationship  between  integrated  mass  and  time  is  more  linear 
than  exponential.  The  interesting  aspect  of  this  plot  is  how  little  mass  has  actually  evaporated.  Finally,  the 
effect  of  vaporization  on  the  shock  wave  Mach  number  is  shown  in  Figure  17. 

The  small  amount  of  vaporization  seen  in  Fig.  16  is  due  primarily  to  the  low  temperatm  e  of  the  shocked 
driven  gases  (about  400  K).  The  vaporization  reduces  this  temperature  quite  rapidly  to  the  limiting  temper¬ 
ature  for  vaporization  (for  this  study,  the  limiting  temperature  is  set  to  310K).  This  is  a  limitation  of  the 
present  simulations,  and  only  studies  with  stronger  shocks  will  show  the  effect  of  appreciable  vaporization. 

Also  of  interest  is  the  considerably  lower  amount  of  vaporization  for  the  equilibrium  model  as  opposed 
to  the  empirical-/S  model.  There  are  several  reasons  for  this  result,  some  of  which  are  physical  and  others 
are  limitations  of  the  model.  First,  the  amount  of  vaporization  for  the  equilibrium  model  is  closely  coupled 
to  the  temperature  of  the  droplet.  At  low  temperatures,  the  amount  of  vaporization  is  minimal,  and  much 
of  the  heat  flux  to  a  droplet  is  used  to  increase  the  temperature  of  the  droplet.  Only  after  the  droplet  is 
heated  will  the  vaporization  rate  become  large.  This  is  physically  realistic  and  is  an  important  improvement 
over  the  empuical-^S  model.  However,  because  we  are  using  the  infinite  conductivity  model,  we  once  again 
artificially  limit  the  surface  temperature  of  the  droplet  by  requiring  the  entire  droplet  to  be  at  a  constant 
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(C)  (d) 

Figure  14:  Snapshot  of  density  solutions  for  the  representative  water  droplet  case  at  1  ms  intervals  for  gas- 
phase.  Mo  =  1.49,  r)  =  0.63.  Water  droplets  initially  evenly  dispersed  from  25  /im  to  30  /lim.  No  (solid), 
temperature-dependent  (dashed),  and  equilibrium  (dash-dot)  vaporization. 
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Figure  15:  Snapshot  of  solutions  for  the  representative  water  droplet  case  at  1  ms  intervals  for  gas-phase. 
Mo  =  1.49,  7}  =  0.63.  Water  droplets  initially  evenly  dispersed  from  25  /im  to  30  ^m.  No  (solid), 
temperature-dependent  (dashed),  and  equilibrium  (dash-dot)  vaporization. 
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Time,  s 

Figure  16:  Total  integrated  mass  of  water  in  individual  sections  for  the  representative  water  droplet  cases. 
Mo  =  1.49, 1)  =  0.63.  Water  droplets  initially  evenly  dispersed  from  25  /xm  to  30  /im. 

temperature.  Near  the  shock  front  (where  vaporization  will  be  most  important),  however,  the  outer  region  of 
the  droplet  will  be  at  a  higher  temperature  and  allow  for  more  vaporization  than  the  interior  regions  of  the 
droplet,  and  thus  the  vaporization  rate  will  be  higher  than  that  predicted  by  the  infinite  conductivity  model. 
Also,  because  of  the  extremely  high  velocities  encountered  within  the  shock,  the  correlations  currently  used 
for  convective  effects  may  be  inadequate.  Therefore  we  feel  that  the  equilibrium  model  represents  a  lower 
bound  on  vaporization,  while  the  empirical-^S  model  serves  as  a  useful  upper  bound. 

6.4  Water  droplet  size  effects 

Two  parametric  studies  were  done  to  examine  the  effect  of  droplet  size  on  shock  mitigation.  The  first 
study  keeps  the  mass-loading  of  water  droplets  constant,  while  varying  the  droplet  size.  This  has  the  effect  of 
also  changing  the  number  density  of  the  droplets.  The  second  set  of  cases  keeps  the  number  density  constant 
while  varying  the  droplet  size.  The  results  (in  terms  of  shock  wave  Mach  number)  are  shown  in  Figure  18. 
As  can  be  seen  comparing  the  two  figures,  die  equilibrium  value  for  the  shock  wave  Mach  number  appears 
to  be  closely  tied  to  the  mass-loading  and  relatively  independent  of  droplet  size  (apparent  when  comparing 
the  12.5-15  fim  and  the  25-30  /rm  results),  even  though  close  to  the  diaphragm  the  mitigation  curves  for  the 
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Figure  17:  Effect  of  vaporization  on  shock  wave  Mach  number  for  the  representative  water  droplet  cases. 
Mo  =  1.49,  T]  =  0.63.  Water  droplets  initially  evenly  dispersed  from  25  pm  to  30  pm. 

constant  number-density  cases  appear  to  more  closely  resemble  each  other.  This  may  imply  that  the  initial 
region  may  be  more  dependent  on  vaporization  processes,  whereas  the  equilibrium  region  is  more  dependent 
on  the  mass-loading,  however,  the  results  are  not  conclusive  in  this  respect. 

7  Shock  Wave  Decay  in  Seeded  Gases 

One  of  the  key  observations  made  in  the  section  above  is  that  maximum  overpressure  is  greater  for  the 
“mitigated”  cases  than  the  unmitigated  cases.  Therefore,  one  may  think  that  injecting  glass  particles  into 
the  flow  is  rather  poor  at  mitigating  the  effects  of  the  blast  overpressure.  However,  we  now  examine  some 
cases  with  much  smaller  driver  sections.  These  cases  have  a  much  shorter  driver  section  of  0.5  m,  while  the 
seeded  driven  section  and  the  unseeded  driven  section  remain  the  same  at  3.4  m  and  0.1  m  respectively.  In 
these  cases,  the  reflected  rarefaction  waves  play  an  important  role  in  the  shock  wave  dynamics. 

7.1  Representative  case  with  glass  particles 

A  representative  solution  was  first  run  with  identical  gas  and  particle  conditions  as  the  representative 
solution  with  the  long  driver  section.  The  seeded  region  contains  a  nearly  mono-dispersed  flow  with  one 
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Figure  18:  Effect  of  droplet  size  on  shock  wave  Mach  number.  Mass-loading  of  water  droplets  is  held 
constant  at  0.63  (a),  number-density  is  held  constant  G>).  Mg  =  1.49.  Empirical-^S  vaporization  model  is 
used. 


Table  2:  Simulations  computed  with  the  short  driver  section  geometry.  Length  of  driver  section  is  0.5  m  for 
all  simulations. _ 


Section  Length 

Mass 

Number 

Particle 

Vaporization 

Case  # 

Dispersed-Phase 

Unseeded 

Seeded 

Loading 

Sections 

Size  (/xm) 

Model 

1 

None 

0.1 

3 

0 

NA 

NA 

2 

Glass 

0.1 

3 

0.63 

1 

26-28 

3 

None 

0.5 

10 

0 

NA 

NA 

4 

Glass 

0.5 

10 

0.63 

NA 

0 

5 

Glass 

0.5 

10 

0.63 

1 

4-6 

6 

Glass 

0.5 

10 

0.63 

1 

9-11 

7 

Glass 

0.5 

10 

0.63 

1 

19-21 

8 

Glass 

0.5 

10 

0.63 

1 

39-41 

9 

Glass 

0.5 

10 

0.63 

1 

59-61 

10 

Glass 

0.5 

10 

0.63 

1 

99-101 

11 

Water 

0.1 

3 

0.63 

6 

0-30 

None 

12 

Water 

0.1 

3 

0.63 

6 

0-30 

Empirical-;0 

13 

Water 

0.1 

3 

0.63 

6 

0-30 

Equilibrium 

14 

Water 

0.1 

3 

0.63 

6 

0-60 

Empirical-;0 

15 

Water 

0.1 

3 

0.63 

6 

0-15 

Empirical-j8 

16 

Water 

0.1 

3 

0.63 

6 

0-120 

Empirical-;0 

17 

Water 

0.1 

3 

0.63 

12 

0-30 

Empirical-yS 

18 

Water 

0.1 

3 

1.26 

6 

0-51.3 

Empirical-;S 

19 

Water 

0.1 

3 

0.315 

6 

0-17.54 

Empirical-|0 
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section  and  a  diameter  range  from  26  ;xm  to  28  /lim.  The  mass-loading  for  all  of  these  cases  is  Jp  — 

0.63.  As  before,  we  show  the  gas-phase  solution  profiles  at  1  ms  intervals  for  the  mitigated  and  unmitigated 
cases  in  Figure  19.  For  the  unmitigated  simulations,  we  see  the  classical  shock-tube  behavior.  Even  before 
1  ms,  the  rarefaction  wave  has  reflected  off  the  back  wall.  The  rarefaction  wave  then  tends  to  reduce  the 
pressure  just  behind  the  shock  wave  until  it  finally  reaches  the  shock  front  and  affects  the  shock  wave.  The 
rarefaction  wave  also  lowers  the  temperature  behind  the  shock  wave  and  substantially  slows  down  the  driver 
gases. 

Comparing  the  unmitigated  (dashed)  and  the  mitigated  (solid)  pressure  results  in  Figure  19c,  the  overall 
pressure  is  clearly  reduced  for  the  mitigation  simulations  away  from  the  initial  diaphragm  (after  approxi¬ 
mately  X  =  1.2  m).  Because  of  the  shape  of  the  pressure  profiles  for  the  mitigated  studies,  the  rarefaction 
wave  has  less  far  to  travel  to  relieve  the  maximum  overpressure  for  the  shock  front,  and  thus  reduces  the 
maximum  overpressure  much  more  quickly  than  the  unmitigated  case.  The  same  is  true  for  the  momentum, 
density,  and  velocity  within  the  shock  front. 

Figure  20  shows  the  gas-phase  and  particle-phase  profiles  for  the  representative  solution.  The  particle- 
phase  profiles  tend  to  follow  the  same  behavior  as  for  the  mitigation  cases  with  long  driver  sections.  One 
difference  is  found  in  the  temperature  solutions.  Figure  20b.  For  the  long  driver  sections,  the  particles  reach 
a  high-temperature  only  within  the  driver  gases.  For  the  mitigated  cases,  the  particle  temperature  increases 
past  the  contact  discontinuity,  but  then  this  high  temperature  moves  into  the  region  of  flow  with  the  driven 
gas.  The  reason  for  this  is  that  the  driver  gas  is  being  slowed  by  the  rarefaction  wave.  Because  the  particles 
adjust  to  the  gas  velocity  slowly  (as  shown  in  Figure  20d),  the  particles  move  from  the  driver  gas  back  into 
the  driven  gas.  This  is  also  why  we  see  a  higher  velocity  for  the  mitigated  cases  away  from  the  initial  shock, 
such  as  in  Figure  19f  for  the  4  ms  (red)  results. 

Close  to  the  diaphragm,  before  the  rarefaction  wave  affects  the  shock  front  (at  a  time  of  1  ms),  the 
pressure  for  the  mitigated  cases  is  slightly  greater  than  the  unmitigated  cases.  Further  downstream,  how¬ 
ever,  shown  in  Figure  19c,  the  overpressure  does  decrease  for  the  mitigation  cases.  To  quantify  better  the 
resulting  overpressure  decrease  due  to  the  glass  particles  for  specific  distances  from  the  diaphragm,  pressure 
sensors  were  placed  at  5  cm  intervals  in  the  seeded  driven  section  of  the  shock  tube.  Taking  the  maximntn 
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Figure  19:  Snapshot  of  representative  solution  at  1  ms  intervals  for  unmitigated  (solid)  and  mitigated 
(dashed)  shocks.  Mq  =  1.49,  r}  =  0.63, 26-28  //m  particles  with  small  driver  section  geometry. 
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Figure  20:  Snapshot  of  representative  solution  at  1  ms  intervals  for  gas-phase  (solid)  and  dispersed-phase 
(dashed)  for  the  mitigating  case.  Mo  =  1.49,  t]  =  0.63,  26-28  fim  particles  with  small  driver  section 
geometiy. 
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Figure  21:  Maximum  overpressure  for  representative  solution  as  a  function  of  sensor  location.  Mo  =  1.49, 
rj  =  0.63,  26-28  pm  particles  with  small  driver  section  geometry. 


overpressure  at  each  point  and  plotting  it  versus  sensor  location,  the  amount  of  pressure  mitigation  observed 
for  a  given  distance  from  the  initial  diaphragm  is  quantified.  This  is  plotted  in  Figure  21  to  give  a  good  idea 
of  the  amount  of  mitigation  to  be  expected  for  the  given  geometry. 

7.2  Glass  particle  size  effects 

Important  to  this  study  is  a  determination  of  the  effect  of  particle  size  on  the  mitigation  characteristics 
of  shock  waves  in  seeded  gases.  A  parametric  study  of  particle  size  was  done  using  the  same  set  up  as 
described  above  for  the  representative  solution.  The  size  of  the  driven  section  is  extended  to  10  meters  such 
that  we  capture  a  larger  portion  of  the  decay  region  of  the  shock  wave.  The  diaphragm  is  located  at  axial 
location  x  =  -0.5  m  such  that  the  unseeded/seeded  interface  is  now  at  the  origin. 

The  range  of  sizes  considered  for  the  parametric  study  varied  from  infinitely  small  particles  (equivalent- 
gas  formulation)  to  100  pm  particles.  The  results  are  divided  into  two  groups:  small  particles  (particles 
smaller  than  around  40  pm)  and  large  particles  (particles  larger  than  around  40  pm).  The  pressure  profile  at 
five  ms  intervals  is  shown  in  Figure  22  for  several  different  particle  sizes.  Fig.  22a  shows  the  smaller  particle 
group,  and  Fig.  22b  shows  the  larger  particle  group.  With  the  smaller  driver  section,  attenuation  starts  very 
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Figure  22:  Snapshots  of  pressure  profiles  at  5  ms  intervals  for  particle  size  study.  Mo  =  1.49,  t)  =  0.63, 
with  the  small  driver  section  geometry  (0.5m)  and  an  extended  seeded  driven  section  (10m). 

soon  after  the  diaphragm  breaks.  The  dotted  line  for  both  plots  indicate  an  unmitigated  case.  Attenuation  of 
the  unmitigated  shock  wave  begins  about  2  m  downstream  from  the  seededAmseeded  interface!  (about  2.5  m 
from  the  diaphragm).  Attenuation  of  the  mitigated  shocks  begins  before  that,  as  we  saw  in  flie  representative 
case.  Looking  at  the  smaller  particles,  as  the  size  increases,  the  peak  pressure  that  is  sharp  for  the  infinitely 
small  particles  becomes  rounded  and  more  diffuse  for  the  larger  particles.  This  is  especially  apparent  for  the 
largest  particles  shown  in  Fig.  22a,  which  shows  considerably  more  overpressure  reduction  than  the  smaller 
particle  cases,  but  less  reduction  in  the  shock  wave  velocity  . 

For  the  larger  droplets,  a  different  trend  is  evident.  Early  in  the  calculation  (at  0.01  seconds)  the  pressure 
profiles  for  the  largest  particles  (60  fjtm  and  100  /rm)  still  resembles  the  unmitigated  pressure  wave  (although 
the  maximum  pressure  is  reduced  substantially).  At  this  point,  the  flow  is  still  transitioning  from  the  un¬ 
mitigated  profile  to  the  mitigated  equilibrium  profiles.  During  this  transition,  the  pressure  is  not  reduced  as 
fast  as  the  smaller  particle  cases.  However,  once  the  transition  is  complete,  the  pressure  reduction  becomes 
much  greater,  and  is  more  similar  to  the  smaller  particle  cases. 

The  above  2  sets  of  results  suggest  diat  there  is  an  optimum  droplet  size.  Another  way  of  look  at  the 
results  is  given  in  Figure  23.  These  plots  show  the  maximum  overpressure  felt  at  any  point  downstream  of 
the  seeded/unseeded  section  interface.  Again,  the  plot  is  divided  into  two  groups  dependent  on  the  particle 
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Figure  23:  Maximum  overpressure  as  a  function  of  sensor  location.  Mo  =  1.49,  r)  =  0.63,  with  the  small 
driver  section  geometry  (0.5m)  and  an  extended  seeded  driven  section  (10m). 

size  for  clarity  of  presentation.  The  graph  of  smaller  droplets  (Figure  23a),  shows  that  increasing  the  size 
clearly  increases  the  amount  of  pressure  reduction  up  to  about  40  urn.  At  this  point,  the  plot  becomes  much 
more  difficult  to  interpret  because  during  the  initial  phases  of  decay,  the  shock  wave  is  still  transitioning  to 
its  equilibrium  setup  with  respect  to  the  particles.  Thus  the  deceleration  of  the  shock  wave  is  slower,  until  it 
completes  the  transition  at  which  point  it  increases.  This  effect  is  more  clearly  shown  for  the  larger  particles 
in  Fig.  23b. 

7.3  Representative  case  with  water  droplets 

As  before,  we  first  examine  the  effect  of  the  vaporization  models  on  the  decay  characteristics  of  the 
shock  wave  for  a  representative  case.  This  representative  solution  uses  the  same  conditions  as  the  previous 
representative  solution  for  glass  particles,  except  that  water  droplets  are  used  instead  of  glass  particles.  The 
density,  pressure,  and  temperature  profiles  are  shown  in  Figure  24.  As  with  the  glass  particles,  downstream 
from  the  diaphragm  the  pressure  is  reduced  for  mitigated  cases  compared  with  the  unmitigated  case,  with 
the  greatest  reductions  in  maximum  overpressure  seen  with  vaporization.  Also  interesting  to  note  is  that 
the  temperature  goes  below  300  K  for  both  vaporization  cases.  This  is  because  the  droplet  vaporization 
tends  to  reduce  the  temperature  drastically  before  the  rarefaction  wave  reflects  and  penetrates  into  the  driven 
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gases,  as  demonstrated  in  Figure  15b  for  steady  shock  waves.  When  the  rarefaction  wave  penetrates  into 
the  driven  gases,  the  temperature  is  further  reduced,  at  which  point  it  goes  below  300  K.  Currently  there  is 
no  condensation  model  in  the  program,  but  clearly  some  condensation  may  occur  at  this  point  raising  the 
temperature.  The  total  integrated  mass  in  several  sections  is  shown  in  Figure  25. 

In  terms  of  mitigation,  we  again  have  pressure  sensors  spaced  25  cm  apart.  Plotting  the  maximum 
overpressure  for  each  sensor  in  Figure  26  and  comparing  the  non-mitigating  and  mitigating  cases,  we  clearly 
see  the  effect  of  vaporization  on  the  maximum  overpressure  felt  at  any  location.  Interesting  to  note  is  that 
the  equilibrium  vaporization  model  initially  gives  results  similar  to  the  non-evaporating  results;  however, 
by  the  end  of  the  domain  the  maximum  overpressure  is  close  to  the  empirical-^  vaporization  results.  This 
is  because  the  equilibrium  model  produces  much  smaller  vaporization  rates  initially,  until  the  surface  of 
the  droplets  increase  in  temperature;  however,  as  the  shock  wave  continues  to  travel  the  total  amount  of 
vaporization  tends  to  approach  the  amount  vaporized  with  the  empirical-^S  model  (Figure  25.) 

7.4  Water  droplet  size  effects 

I 

The  next  set  of  simulations  looked  at  the  effect  of  droplet  size  on  mitigation  of  the  maximum  overpres¬ 
sure.  As  with  the  previous  studies  on  droplet  size,  we  do  two  sets  of  simulations  and  compared  the  results. 
All  of  these  cases  used  the  empirical-/?  vaporization  model.  The  first  set  kept  the  mass-loading  constant  at 
0.63  and  varied  the  droplet  size  from  12.5-15  jim  to  100-120  //m.  The  second  set  of  computations  kept  the 
number  density  constant  for  all  cases  and  varied  the  mass-loading  with  the  particle  size. 

The  results  for  maximum  overpressure  as  a  function  of  location  for  both  series  of  computations  are 
shown  in  Figure  27.  Because  the  shock  wave  is  fairly  weak,  the  temperature  of  the  shocked  driven  gases 
only  reaches  about  400K,  and  this  is  rapidly  cooled  by  evaporation  of  only  a  small  part  of  the  water  avail¬ 
able.  Because  of  this,  the  parametric  study  with  the  mass-loading  constant  (Fig.  27a)  follows  closely  the 
characteristics  of  the  glass  particle  study  presented  earlier  in  this  report. 

The  parametric  study  keeping  the  number  density  constant  as  the  droplet  size  and  mass-loading  varies 
is  shown  in  Fig.  27b.  The  idea  behind  this  study  is  that  for  the  smaller  droplets,  a  larger  percentage  of  the 
available  water  is  able  to  vaporize,  thus  the  reduction  in  pressure  (relative  to  the  larger  droplets  and  larger 
mass-loading)  is  more  similar  even  though  much  less  water  is  used.  From  the  results,  it  is  difficult  to  tell 
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Figure  24:  Snapshot  of  solution  for  representative  water  droplet  case  at  1  ms  time  intervals.  Mo  =  1.49, 
T)  =  0.63,  with  short  driver  section  geometry.  Water  droplets  initially  evenly  dispersed  from  25  ^m  to  30 
/xm.  No  (solid),  temperature-dependent  (dashed),  and  equilibrium  (dash-dot)  vaporization. 
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Figure  25:  Total  integrated  mass  of  water  in  individual  sections  for  representative  water  droplet  case.  Mo  = 
1.49,  ri  =  0.63,  with  short  driver  section  geometry.  Water  droplets  initially  evenly  dispersed  from  25  //m  to 
30  fiva.  No  (solid),  temperature-dependent  (dashed),  and  equilibrium  (dash-dot)  vaporization. 
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Figure  26:  Maximum  overpressure  as  a  function  of  sensor  location.  Mo  =  1.49,  rj  =  0.63,  witii  the  small 
driver  section  geometry. 
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(a)  (b) 

Figure  27:  Maximum  overpressure  for  different  droplet  sizes  with  mass-loading  held  constant  at  77  =  0.63 
(left),  and  with  number  density  held  constant  (right).  Mo  =  1.49,  with  the  small  driver  section  geometry. 


whether  there  is  any  benefit  from  using  smaller  droplets  and  a  smaller  mass-loading.  Part  of  the  reason  for 
this  is  that  it  is  not  appropriate  to  keep  the  number  density  constant.  A  more  appropriate  study  would  keep 
the  surface  area  constant  between  different  droplet  sizes,  thus  the  amount  of  vaporization  is  similar.  For 
our  case,  there  will  still  be  more  total  vaporization  for  the  larger  droplet  cases,  even  though  the  percentage 
of  water  evaporated  is  smaller.  Apart  from  this,  these  cases  suffer  from  the  same  problem  as  the  other 
vaporization  studies;  namely,  that  the  shock  wave  is  too  weak  to  create  high  temperatures  in  the  shocked 
gases,  limiting  the  amount  of  vaporization  that  will  occur. 

8  Summary  and  Conclusions 


Because  of  the  number  of  simulations  conducted  for  this  report,  we  first  summarize  our  results  before 
presenting  conclusions  discussing  the  relevance  to  blast  mitigation.  The  report  started  with  a  summary  of 
shock  relations  for  transmission  of  shock  waves  through  particle-laden  flows  for  infinitely  small  particles. 
From  this  analysis,  it  was  found  that  for  most  cases,  transmission  of  the  shock  wave  results  in  compression 
of  the  particle-laden  gases  by  a  reflected  shock  wave,  an  increase  in  the  overpressure  behind  the  shock 
wave  and  a  decrease  in  the  shock  wave  velocity.  The  increase  in  overpressure  is  caused  by  compression  of 
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the  gases  due  to  lower  gas  velocities.  This  analysis,  however,  did  not  account  for  particle  size  effects  or 
vaporization  effects,  which  requires  a  numerical  treatment. 

The  sectional  approach  for  dispersed-phase  numerical  simulations  was  consequently  developed  for  a 
convective  flow-field.  This  approach  allows  the  dispersed-phase  to  be  treated  by  continuum  conservation 
equations.  Cross-coupling  terms  between  the  gas  and  dispersed-phase,  and  between  the  different  sections 
in  the  dispersed-phase  were  developed  using  single-droplet  models.  The  infinite  conductivity  model  was 
used  for  calculating  the  energy  of  the  particles,  and  the  equilibrium  vaporization  model  was  used  for  cal¬ 
culating  vaporization  rates  with  corrections  to  account  for  convective  effects.  The  sectional  approach  was 
then  incorporated  into  a  high-speed  computational  fluid  dynamics  code  based  on  the  flux-corrected-transport 
algorithm. 

The  numerical  methodology  was  first  validated  using  experimental  data.  The  experiments  investigated 
the  transition  of  the  shock  wave  from  unseeded  gases  to  gases  seeded  with  glass  particles.  Good  agreement 
was  found  with  the  results  for  different  Mach  number  shock  waves  and  different  mass-loading  of  particles. 
Interestingly,  it  was  found  that  better  agreement  could  be  found  at  the  beginning  of  the  transition  region 
for  simulations  that  did  not  include  heat  transfer,  although  these  simulations  generally  did  not  get  the  equi¬ 
librium  shock  wave  Mach  number  correct.  This  is  most  likely  caused  by  the  infinite  conductivity  model, 
which  keeps  the  particle  surface  temperature  artificially  low  initially,  thus  overpredicting  heat  transfer  to  the 
particles. 

Further  simulations  were  conducted  to  examine  the  mitigation  characteristics  of  steady  and  decaying 
shock  waves  impinging  on  seeded  flows.  First,  we  simulated  the  effects  of  steady  shock  waves,  produced 
using  shock  tubes  with  long  driver  sections.  We  described  both  the  effect  of  non-vaporizing  particles  and 
vaporizing  droplets  on  the  pressure  and  velocity  of  the  shock  wave.  We  found  that  although  seeding  the 
flow  generally  causes  the  shock  wave  to  slow  down,  the  pressure  developed  behind  the  shock  wave  actually 
increases  in  particle-laden  flows.  For  infinitely  small  particles,  the  square  wave  form  is  maintained,  however, 
for  finite  sized  particles  the  shock  structure  is  spread  out  such  that  the  maximum  overpressure  is  developed 
downstream  of  the  actual  shock  front.  With  vaporization,  the  temperature  is  reduced  considerably.  The 
pressure,  however,  is  reduced  only  slightly  because  the  density  of  the  gas  has  now  also  increased  significantly 
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due  to  the  water  vapor.  For  the  cases  presented  in  this  report,  the  amount  of  vaporization  is  quite  small  due  to 
the  weak  shocks  we  investigated,  which  result  in  only  moderate  temperatures  for  the  shocked  driven  gases. 

Another  set  of  simulations  examined  decaying  shock  waves,  produced  using  shock  tubes  with  short 
driver  sections.  Again,  we  looked  at  the  effect  of  different  particle  sizes  on  the  velocity  of  the  shock  wave  and 
the  pressure  developed  behind  the  shock  front.  Unlike  the  steady  shock  waves,  the  maximum  overpressure 
for  decaying  shock  waves  is  reduced  for  the  particle-seeded  simulations.  This  effect  is  caused  because  the 
droplets  tend  to  disperse  energy  away  from  the  sharp  peak.  As  the  particle  size  increases,  this  effect  becomes 
more  pronounced.  However,  there  is  a  conflicting  effect.  If  the  shock  wave  has  not  fiilly  transitioned  to 
the  final  shock  structure  corresponding  to  the  equilibrium  condition,  the  particles  are  not  as  effective  in 
dispersing  the  energy  away  from  the  sharp  peak.  This  results  in  higher  maximum  overpressures  compared 
to  smaller  droplets.  The  effect  of  vaporization  for  these  cases  is  to  further  reduce  the  maximum  overpressure 
felt  at  specific  locations.  The  degree  to  which  the  maximum  overpressure  is  reduced  depends  on  the  amount 
of  vaporization  that  occurs.  For  the  simulations  done,  we  found  a  small  overall  reduction  in  maximum 
overpressure  and  a  somewhat  larger  reduction  in  the  velocity  of  the  shock  wave.  As  for  the  previous  studies, 
the  small  effect  of  vaporization  on  the  maximum  overpressures  is  due  to  the  low  temperature  of  the  shocked 
driven  gases. 

Additional  studies  are  required  to  further  quantify  the  effects  of  water  droplets  on  blast  mitigation  and 
how  to  effectively  mitigate  them.  In  particular,  vaporization  needs  to  be  more  closely  studied.  These  sim¬ 
ulations  indicated  that  vaporization  decreases  the  air  temperature  significantly,  but  only  has  a  small  effect 
on  the  shock  overpressure.  Because  of  the  low  temperatures  in  the  shocked  driven  gases  for  these  cases, 
vaporization  had  only  a  secondary  role  in  the  mitigation  studies  conducted  above.  A  stronger  shock  wave  is 
essential  to  allow  the  vaporization  to  have  a  more  critical  role  in  mitigating  the  shock  wave  produced  from 
either  a  blast  or  a  shock  tube. 

Blasts  from  explosives  are  in  many  ways  similar  to  the  shock  tube  flows  used  for  this  study.  The  initial 
explosive  is  analogous  to  the  shock  tube  driver  section,  and  acts  to  form  the  shock  wave  after  the  energy 
is  released.  For  blasts,  because  of  the  small  size  of  the  explosive  and  the  fact  that  most  blasts  spherically 
emanate  from  the  initial  explosive  (unlike  shock  tubes,  which  are  usually  planar),  the  shock  wave  for  the 
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blast  almost  instantly  starts  decaying.  From  the  simulations  with  the  planar  shock  wave,  we  would  expect  a 
particle-laden  flow  to  effectively  mitigate  this  type  of  shock  wave  since  it  is  similar  to  the  decaying  shock 
waves  discussed  in  this  report. 

Although  shock  tubes  and  explosive  blasts  share  many  similarities,  they  also  have  significant  differences. 
In  addition  to  the  previously  mentioned  difference  between  planar  and  spherical  shock  wave,  explosive  blasts 
have  density  ratios  that  are  significantly  greater  (by  orders  of  magnitude)  than  these  produced  in  a  shock 
tube.  Also,  many  explosives  are  oxygen-deficient,  and  thus  have  reactions  occurring  within  the  driver  gases. 
Penetration  of  the  droplets  into  the  driver  gas  region  will  lead  to  additional  mitigating  effects.  Although 
numerical  simulations  may  guide  experimental  work  for  determining  mitigation  of  explosive  blasts  from 
water  mist,  experimental  work  is  crucial  in  determining  the  validity  of  the  simulation  technique  and  for 
refining  the  models  to  the  extreme  environments  found  for  explosive  blasts. 
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